{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "c3062d31",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from functools import reduce"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "fc19d05b",
   "metadata": {},
   "outputs": [],
   "source": [
    "a=np.array([3,4,5,6,7])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "a6f6d64b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0, 1, 2, 4]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a=[i for i in range(5) if i!=3]\n",
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "30c5ad6a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ans=reduce(lambda x,y: x*y, [i for i in range(5) if i!=3])\n",
    "ans"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5096eb94",
   "metadata": {},
   "source": [
    "## Lagrange Polynomial Interpolation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 229,
   "id": "5612caa9",
   "metadata": {},
   "outputs": [],
   "source": [
    "class Lagrange_Polynomial_Interpolation(object):\n",
    "    def __init__(self, x, y):\n",
    "        self.N=x.shape[0]\n",
    "        self.lx=[]\n",
    "        self.x=np.copy(x)\n",
    "        self.y=np.copy(y)\n",
    "    \n",
    "    def calculate_lk(self, x, k,):\n",
    "        ta=reduce(lambda x,y: x*y, [x-self.x[i] for i in range(self.N) if i!=k])\n",
    "        tb=reduce(lambda x,y: x*y, [self.x[k]-self.x[i] for i in range(self.N) if i!=k])\n",
    "        return ta/tb\n",
    "    \n",
    "    def predict(self, x):\n",
    "        return np.sum([self.y[k]*self.calculate_lk(x,k) for k in range(self.N)])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "id": "091896c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "x=np.array([0.32,0.34,0.36])\n",
    "y=np.array([0.314567,0.333487, 0.352274])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 879,
   "id": "e05a5145",
   "metadata": {},
   "outputs": [],
   "source": [
    "LPI=Lagrange_Polynomial_Interpolation(x,y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 880,
   "id": "db86bbf6",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.3303743620375"
      ]
     },
     "execution_count": 880,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x0=0.3367\n",
    "LPI.predict(x0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 881,
   "id": "83047069",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-1.7048187195278786e-07"
      ]
     },
     "execution_count": 881,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sin(x0)-LPI.predict(x0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "90b5b9c5",
   "metadata": {},
   "outputs": [],
   "source": [
    "F=lambda x: np.sin(x**2)*np.exp(-2*x)+x**2+4*x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "946c0667",
   "metadata": {},
   "outputs": [],
   "source": [
    "x=np.arange(0,10,1)\n",
    "y=F(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "id": "2cf5d064",
   "metadata": {},
   "outputs": [],
   "source": [
    "LPI=Lagrange_Polynomial_Interpolation(x,y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "id": "cad0129e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "id": "11dd1d5f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  0.        ,   5.11388071,  11.98613868,  21.00102154,\n",
       "        31.99990342,  44.99999399,  59.99999391,  76.99999921,\n",
       "        96.0000001 , 116.99999999])"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "id": "f3988001",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "65.49748708199867"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "LPI.predict(6.3366)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "id": "4cdf4ed1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "65.49890155056632"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F(6.3366)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 194,
   "id": "6c6176a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "class Successive_Linear_Interpolation(object):\n",
    "    def __init__(self, x, y, base_order):\n",
    "        self.N=x.shape[0]\n",
    "        if base_order is not None:\n",
    "            assert base_order<=self.N,'The base_order should less than the number of points'\n",
    "        self.lx=[]\n",
    "        self.x=np.copy(x)\n",
    "        self.y=np.copy(y)\n",
    "        self.base_order=base_order if base_order is not None else self.N\n",
    "        self.I_set=[[] for k in range(self.base_order)]\n",
    "        self.calculated=False\n",
    "        \n",
    "    def calculate_I(self, x, eps=1e-6):\n",
    "        value=[[] for k in range(self.base_order)]\n",
    "        # I_k=lambda k,j, x: self.y[k]+(self.y[j]-self.y[k])/(self.x[j]-self.x[k])*(x-self.x[k])\n",
    "        F=lambda i,j,m,n,s,r, x: self.I_set[i][j]+(self.I_set[m][n]-self.I_set[i][j])/(self.x[s]-self.x[r])*(x-self.x[r])\n",
    "        if not self.calculated:\n",
    "            self.I_set[0].append(self.y[0])\n",
    "\n",
    "            for i in range(1,self.base_order):\n",
    "                self.I_set[i].append(self.y[i])\n",
    "                self.I_set[i].append(F(i-1,0,i,0,i,i-1, x))\n",
    "                for k in range(len(self.I_set[i-1])-1):\n",
    "                    self.I_set[i].append(F(i-1,1+k, i, 1+k, i, i-2-k, x))\n",
    "\n",
    "                if i>1:\n",
    "                    if abs(last-self.I_set[i][-1])<=eps:\n",
    "                        return self.I_set[i][-1]\n",
    "                    else:\n",
    "                        last=self.I_set[i][-1]\n",
    "                else:\n",
    "                    last=self.I_set[i][-1]\n",
    "        \n",
    "        \n",
    "        \n",
    "        return self.I_set[-1][-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "91ea3190",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.  0.2 0.3 0.5 0.6] [0.      0.20134 0.30452 0.5211  0.63665]\n"
     ]
    }
   ],
   "source": [
    "x=np.array([0,0.2,0.3,0.5,0.6])\n",
    "y=np.array([0,0.20134,0.30452,0.52110,0.63665])\n",
    "print(x,y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 203,
   "id": "ad9692a1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.32 0.34 0.36]\n",
      "[0.314567 0.333487 0.352274]\n"
     ]
    }
   ],
   "source": [
    "print(x)\n",
    "print(y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 207,
   "id": "e059d88c",
   "metadata": {},
   "outputs": [],
   "source": [
    "SLI=Successive_Linear_Interpolation(x,y,None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 208,
   "id": "a935af66",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.3303743620375"
      ]
     },
     "execution_count": 208,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SLI.calculate_I(0.3367)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 167,
   "id": "9c571b41",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.232294"
      ]
     },
     "execution_count": 167,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y[1]+(y[2]-y[1])/(x[2]-x[1])*(0.23-x[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 829,
   "id": "1a90b327",
   "metadata": {},
   "outputs": [],
   "source": [
    "class mini_I(object):\n",
    "    def __init__(self, x0=0,x1=0,y0=0,y1=0, first=False):\n",
    "        self.x0=x0\n",
    "        self.x1=x1\n",
    "        self.y0=y0\n",
    "        self.y1=y1\n",
    "        self.first=first\n",
    "        self.last_answer={}\n",
    "    def __call__(self,x):\n",
    "        # 递归\n",
    "        if x not in self.last_answer.keys():\n",
    "            answer=self.y0(x)+(self.y1(x)-self.y0(x))/(self.x1-self.x0)*(x-self.x0) if not self.first else self.y0\n",
    "            self.last_answer[x]=answer\n",
    "            return answer\n",
    "        else:\n",
    "            if x in self.last_answer.keys():\n",
    "                return self.last_answer[x]\n",
    "            \n",
    "\n",
    "class Successive_Linear_Interpolation(object):\n",
    "    def __init__(self, x, y, base_order):\n",
    "        self.N=x.shape[0]\n",
    "        if base_order is not None:\n",
    "            assert base_order<=self.N,'The base_order should be less than the number of points'\n",
    "        self.lx=[]\n",
    "        self.x=np.copy(x)\n",
    "        self.y=np.copy(y)\n",
    "        self.base_order=0\n",
    "        self.calculated=False\n",
    "        \n",
    "        # I_k=lambda k,j, x: self.y[k]+(self.y[j]-self.y[k])/(self.x[j]-self.x[k])*(x-self.x[k])\n",
    "        #F=lambda i,j,m,n,s,r, x: self.I_set[i][j](x)+(self.I_set[m][n](x)-self.I_set[i][j](x))/(self.x[s]-self.x[r])*(x-self.x[r])\n",
    "        #F=lambda i,j,x: self.I_set[i-1][j-1](x)+(self.I_set[i][j-1](x)-self.I_set[i-1][j-1](x))/(self.x[i]-self.x[i-j])*(x-self.x[i-j])\n",
    "        \n",
    "        self.F=lambda i,j: mini_I(self.x[i-j],self.x[i],self.I_set[i-1][j-1],self.I_set[i][j-1], first=False)\n",
    "        \n",
    "        \n",
    "    def calculate_I(self, x, eps=1e-6, order=None):\n",
    "        assert order<self.N,'The order should be less than the number of points'\n",
    "        if not self.calculated:\n",
    "            \n",
    "            self.base_order=self.N if order is None else order\n",
    "        \n",
    "            self.I_set=[[]]\n",
    "        \n",
    "            self.I_set[0].append(mini_I(y0=self.y[0], first=True))\n",
    "            \n",
    "            last=0\n",
    "            \n",
    "            for i in range(1,self.base_order+1):\n",
    "                self.I_set+=[[]]\n",
    "                self.I_set[i].append(mini_I(y0=self.y[i], first=True))\n",
    "                for k in range(1,i+1):\n",
    "                    self.I_set[i].append(self.F(i,k))\n",
    "                    \n",
    "                if i>1:\n",
    "                    new=self.I_set[i][-1](x)\n",
    "                    print(i,abs(last-new)<=eps)\n",
    "                    if abs(last-new)<=eps:\n",
    "                        if i<(self.base_order+1)-1:\n",
    "                            print(\"Meet with absotion condition of eps={} and the current order of interpolation is {}\".format(eps, i))\n",
    "                            self.base_order=i\n",
    "#                             self.I_set.pop(-1)\n",
    "                            \n",
    "                        self.calculated=True\n",
    "                        return new\n",
    "\n",
    "                last=self.I_set[i][-1](x)\n",
    "                \n",
    "                print(i,)\n",
    "                self.print(x)\n",
    "                \n",
    "            self.calculated=True\n",
    "        \n",
    "        else:\n",
    "            \n",
    "            if self.base_order>=order:\n",
    "                last=self.I_set[0][0](x)\n",
    "                for i in range(1,order+1):\n",
    "                    new=self.I_set[i][-1](x)\n",
    "                    if abs(last-new)<=eps:\n",
    "                        print(\"Meet with absotion condition of eps={} and the current order of interpolation is {}\".format(eps, i))\n",
    "                        return new\n",
    "                    last=new\n",
    "                return self.I_set[order][order](x)\n",
    "            elif self.base_order<order and order <=self.N:\n",
    "#                 self.I_set+=[[] for _ in range((order+1)-self.base_order)]\n",
    "                last=self.I_set[0][0](x)\n",
    "                for i in range(1,self.base_order+1):\n",
    "                    new=self.I_set[i][-1](x)\n",
    "                    if abs(last-new)<=eps:\n",
    "                        print(\"Meet with absotion condition of eps={} and the current order of interpolation is {}\".format(eps, i))\n",
    "                        return\n",
    "                    last=new\n",
    "                \n",
    "                for i in range(self.base_order+1,(order+1)):\n",
    "                    \n",
    "                    self.I_set+=[[]]\n",
    "                    self.I_set[i].append(mini_I(y0=self.y[i], first=True))\n",
    "                    for k in range(1,i+1):\n",
    "                        self.I_set[i].append(self.F(i,k))\n",
    "                        \n",
    "                    if i>self.base_order:\n",
    "                        new=self.I_set[i][-1](x)\n",
    "                        if abs(last-new)<=eps:\n",
    "                            if i<(order+1)-1:\n",
    "                                print(\"Meet with absotion condition of eps={} and the current order of interpolation is {}\".format(eps, i))\n",
    "                                self.base_order=i\n",
    "                                \n",
    "                            return new\n",
    "\n",
    "                    last=self.I_set[i][-1](x)\n",
    "                    \n",
    "                    \n",
    "                self.base_order=(order+1)\n",
    "                return self.I_set[(order+1)-1][(order+1)-1](x)\n",
    "            \n",
    "#         for i in range(len(self.I_set)):\n",
    "#             for j in range(len(self.I_set[i])):\n",
    "#                 print(i,j,self.I_set[i][j](x))\n",
    "        return self.I_set[self.base_order][self.base_order](x)\n",
    "\n",
    "    def add_point(self,x,y):\n",
    "        M=x.shape[0]\n",
    "        \n",
    "        if self.calculated:\n",
    "            self.x=np.concatenate((self.x,x),axis=0)\n",
    "            self.y=np.concatenate((self.y,y),axis=0)\n",
    "#             self.I_set+=[[] for _ in range(M)]\n",
    "#             for i in range(self.N,self.N+M):\n",
    "#                 self.I_set[i].append(mini_I(y0=self.y[i], first=True))\n",
    "#                 for k in range(1,i+1):\n",
    "#                     self.I_set[i].append(self.F(i,k))\n",
    "\n",
    "            \n",
    "            \n",
    "        else:\n",
    "            self.x=np.concatenate((self.x,x),axis=0)\n",
    "            self.y=np.concatenate((self.y,y),axis=0)\n",
    "#             self.I_set+=[[] for _ in range(M)]\n",
    "        \n",
    "        self.N=self.N+M\n",
    "            \n",
    "    def print(self,x):\n",
    "        for i in range(len(self.I_set)):\n",
    "            print('I_{}:'.format(i),end='\\t')\n",
    "            for j in range(len(self.I_set[i])):\n",
    "                print(self.I_set[i][j](x), end='\\t')   \n",
    "            print('\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "id": "23750329",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "\n",
    "class mini_I(object):\n",
    "    \"\"\"\n",
    "    class for calculating the interpolation of order k+1 from two component of order k\n",
    "    base funtion: ans=y0+(y1-y0)/(x1-x0)*(x-x0)\n",
    "\n",
    "    Author: Junno\n",
    "    Date: 2022-03-01\n",
    "    \"\"\"\n",
    "\n",
    "    def __init__(self, x0=0, x1=0, y0=0, y1=0, first=False):\n",
    "        \"\"\"__init__ \n",
    "        Args:\n",
    "            x0 (float, optional): [component of the base function]. Defaults to 0.\n",
    "            x1 (float, optional): [component of the base function]. Defaults to 0.\n",
    "            y0 (float, optional): [component of the base function]. Defaults to 0.\n",
    "            y1 (float, optional): [component of the base function]. Defaults to 0.\n",
    "            first (bool, optional): [whether the first element of each order list]. Defaults to False.\n",
    "        \"\"\"\n",
    "        self.x0 = x0\n",
    "        self.x1 = x1\n",
    "        self.y0 = y0\n",
    "        self.y1 = y1\n",
    "        self.first = first\n",
    "        # record answer for fast calculation\n",
    "        self.last_answer = {}\n",
    "\n",
    "    def __call__(self, x=0.):\n",
    "        \"\"\"__call__ \n",
    "\n",
    "        Args:\n",
    "            x ([float]): [interpolation node]\n",
    "\n",
    "        Returns:\n",
    "            [answer]: [recursive call or the recorded values]\n",
    "\n",
    "        \"\"\"\n",
    "        if x not in self.last_answer.keys():\n",
    "            # here x0,x1,y0,y1 can be a mini_I class to call for input x\n",
    "            answer = self.y0(x)+(self.y1(x)-self.y0(x))/(self.x1 -\n",
    "                                                         self.x0)*(x-self.x0) if not self.first else self.y0\n",
    "            # record\n",
    "            self.last_answer[x] = answer\n",
    "            return answer\n",
    "        else:\n",
    "            if x in self.last_answer.keys():\n",
    "                return self.last_answer[x]\n",
    "\n",
    "\n",
    "class Successive_Linear_Interpolation(object):\n",
    "    \"\"\"Successive_Linear_Interpolation_2D 逐次线性插值\n",
    "\n",
    "    Author: Junno\n",
    "    Date: 2022-03-01\n",
    "    \"\"\"\n",
    "\n",
    "    def __init__(self, x, y, base_order=0):\n",
    "        \"\"\"__init__ \n",
    "\n",
    "        Args:\n",
    "            x ([np.array]): [x of points]\n",
    "            y ([np.array]): [y of points]\n",
    "            base_order ([int]): [interpolation order]\n",
    "        \"\"\"\n",
    "        self.N = x.shape[0]\n",
    "        assert base_order < self.N, 'The base_order should be less than the number of points'\n",
    "        self.x = np.copy(x)\n",
    "        self.y = np.copy(y)\n",
    "        self.base_order = 0\n",
    "        # flag of calculated\n",
    "        self.calculated = False\n",
    "\n",
    "        # ref:《数值分析》 李庆扬、王能超、易大义（华中科大出版社）\n",
    "        self.F = lambda i, j: mini_I(\n",
    "            self.x[i-j], self.x[i], self.I_set[i-1][j-1], self.I_set[i][j-1], first=False)\n",
    "\n",
    "    def calculate_I(self, x, eps=1e-6, order=1, check=True):\n",
    "        \"\"\"calculate_I calculate y for a given x and order with threshold eps\n",
    "\n",
    "        Args:\n",
    "            x ([float]): [x]\n",
    "            eps ([float], optional): [threshold to stop process]. Defaults to 1e-6.\n",
    "            order ([int], optional): [interpolation order]. Defaults to 1.\n",
    "        \"\"\"\n",
    "        assert order < self.N, 'The order should be less than the number of points'\n",
    "\n",
    "        if not self.calculated:\n",
    "\n",
    "            self.base_order = self.N if order is None else order\n",
    "\n",
    "            self.I_set = [[]]\n",
    "\n",
    "            self.I_set[0].append(mini_I(y0=self.y[0], first=True))\n",
    "\n",
    "            last = 0\n",
    "\n",
    "            for i in range(1, self.base_order+1):\n",
    "                self.I_set += [[]]\n",
    "                self.I_set[i].append(mini_I(y0=self.y[i], first=True))\n",
    "                for k in range(1, i+1):\n",
    "                    self.I_set[i].append(self.F(i, k))\n",
    "\n",
    "                # early stop when the difference between of two results from k order and k+1 order is less than eps\n",
    "                if i > 1:\n",
    "                    new = self.I_set[i][-1](x)\n",
    "                    # print(i, abs(last-new) <= eps)\n",
    "                    if abs(last-new) <= eps:\n",
    "                        if i < (self.base_order+1)-1:\n",
    "                            print(\n",
    "                                \"Meet with absotion condition of eps={} and the current order of interpolation is {}\".format(eps, i))\n",
    "                            self.base_order = i\n",
    "#                             self.I_set.pop(-1)\n",
    "\n",
    "                        self.calculated = True\n",
    "                        return new\n",
    "\n",
    "                last = self.I_set[i][-1](x)\n",
    "\n",
    "                # print(i,)\n",
    "                # self.print(x)\n",
    "\n",
    "            self.calculated = True\n",
    "            \n",
    "            return self.I_set[self.base_order][self.base_order](x)\n",
    "\n",
    "        else:\n",
    "            last = self.I_set[0][0](x)\n",
    "            for i in range(1, order+1):\n",
    "\n",
    "                if i > self.base_order:\n",
    "                    self.I_set += [[]]\n",
    "                    self.I_set[i].append(mini_I(y0=self.y[i], first=True))\n",
    "                    for k in range(1, i+1):\n",
    "                        self.I_set[i].append(self.F(i, k))\n",
    "\n",
    "                new = self.I_set[i][-1](x)\n",
    "                if abs(last-new) <= eps:\n",
    "                    if check:\n",
    "                        print(\n",
    "                            \"Meet with absotion condition of eps={} and the current order of interpolation is {}\".format(eps, i))\n",
    "                    return new\n",
    "                last = new\n",
    "\n",
    "            return self.I_set[order][order](x)\n",
    "\n",
    "    def add_point(self, x, y):\n",
    "        \"\"\"add_point \n",
    "\n",
    "        Args:\n",
    "            x ([np.array]): [x of points]\n",
    "            y ([np.array]): [y of points]\n",
    "        \"\"\"\n",
    "        M = x.shape[0]\n",
    "        self.x = np.concatenate((self.x, x), axis=0)\n",
    "        self.y = np.concatenate((self.y, y), axis=0)\n",
    "        self.N = self.N+M\n",
    "\n",
    "    def print(self, x):\n",
    "        \"\"\"print show I_set for a given x\n",
    "\n",
    "        Args:\n",
    "            x ([float]): [x]\n",
    "        \"\"\"\n",
    "        for i in range(len(self.I_set)):\n",
    "            print('I_{}:'.format(i), end='\\t')\n",
    "            for j in range(len(self.I_set[i])):\n",
    "                print(self.I_set[i][j](x), end='\\t')\n",
    "            print('\\n')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "id": "45e74657",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.  0.2 0.3 0.5 0.6] [0.      0.20134 0.30452 0.5211  0.63665]\n"
     ]
    }
   ],
   "source": [
    "x=np.array([0,0.2,0.3,0.5,0.6])\n",
    "y=np.array([0,0.20134,0.30452,0.52110,0.63665])\n",
    "print(x,y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "id": "16793e9e",
   "metadata": {},
   "outputs": [],
   "source": [
    "SLI=Successive_Linear_Interpolation(x,y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "8d72a39f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.34155135758333344"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SLI.calculate_I(0.3355, order=2, check=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "9bbb747f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.014509882820716802"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SLI.calculate_I(0.3367, order=2)-np.sin(0.3345)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "ac1c2d04",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "I_0:\t0.0\t\n",
      "\n",
      "I_1:\t0.20134\t0.33774785\t\n",
      "\n",
      "I_2:\t0.30452\t0.3411489000000001\t0.34155135758333344\t\n",
      "\n",
      "I_3:\t0.5211\t0.34296295000000004\t0.34196824591666675\t0.3418310896550001\t\n",
      "\n",
      "I_4:\t0.63665\t0.3310202499999999\t0.34154973050000004\t0.3418264738192709\t0.3418285086335215\t\n",
      "\n"
     ]
    }
   ],
   "source": [
    "SLI.print(0.3355)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "id": "c791b6bd",
   "metadata": {},
   "outputs": [],
   "source": [
    "SLI.add_point(np.array([0.38+i*0.02 for i in range(5)]),np.array([np.sin(0.38+i*0.02) for i in range(5)]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "id": "aba81fc1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.32, 0.34, 0.36, 0.38, 0.4 , 0.42, 0.44, 0.46])"
      ]
     },
     "execution_count": 73,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SLI.x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "id": "4ea202a7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Meet with absotion condition of eps=1e-06 and the current order of interpolation is 4\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.35881665121993456"
      ]
     },
     "execution_count": 74,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SLI.calculate_I(0.367, order=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "id": "22f2d57d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "I_0:\t0.314567\t\n",
      "\n",
      "I_1:\t0.333487\t0.35902899999999993\t\n",
      "\n",
      "I_2:\t0.352274\t0.35884945\t0.35881802875\t\n",
      "\n",
      "I_3:\t0.3709204694129827\t0.3588002642945439\t0.35881624964881714\t0.3588166351207401\t\n",
      "\n",
      "I_4:\t0.3894183423086505\t0.3588968520307986\t0.3588171671483885\t0.35881666252362426\t0.35881665121993456\t\n",
      "\n"
     ]
    }
   ],
   "source": [
    "SLI.print(0.367)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 874,
   "id": "ab6faf88",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4"
      ]
     },
     "execution_count": 874,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SLI.base_order"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 590,
   "id": "db7ba729",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4.313996905258577e-07"
      ]
     },
     "execution_count": 590,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "0.3875759142025974-0.3875754828029069"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 270,
   "id": "7404f2c3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.32, 0.34, 0.36])"
      ]
     },
     "execution_count": 270,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 271,
   "id": "ea2af085",
   "metadata": {},
   "outputs": [],
   "source": [
    "a=x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 272,
   "id": "1db17600",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.32, 0.34, 0.36, 0.32, 0.34, 0.36])"
      ]
     },
     "execution_count": 272,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.concatenate((a,x), axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 136,
   "id": "1c9c596c",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.patches as mpatches"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "id": "e2d1903f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(100,)"
      ]
     },
     "execution_count": 118,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x=np.arange(0,1,0.01)\n",
    "F=lambda x:np.cos(2*np.pi*x)+np.sin(x)+np.exp(x**2)\n",
    "x.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 204,
   "id": "4127f9f9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.14 0.17 0.32 0.38 0.44 0.52 0.67 0.75 0.88 0.98]\n"
     ]
    }
   ],
   "source": [
    "xs=np.sort(np.random.choice(x, size=(10), replace=False))\n",
    "print(xs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 242,
   "id": "1fc4a253",
   "metadata": {},
   "outputs": [],
   "source": [
    "ys=F(xs)\n",
    "ys_noise=F(xs)+np.random.randn(xs.shape[0],)*0.01"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 243,
   "id": "a315bede",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x210e120e460>"
      ]
     },
     "execution_count": 243,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAT1UlEQVR4nO3df2xd533f8fdHEoNGdjoHErc4piV2g/fLxVx7nOUs3qB5Bha7Ro0B/sMbGwPOAMJWPHhAgXatAAfBQGD7R4gdwxKINE2NEgmKJsi8QN5gIHVjIbMzyr9iR92gdSRN2IBpZbUrqagn6bs/7nVCUZfkpXTJSx6+X8AFz3nOw3O/eiB9ePSc5/CmqpAkbX7b+l2AJKk3DHRJaggDXZIawkCXpIYw0CWpIXb06413795dw8PD/Xp7SdqUjh8//l5VDXY61rdAHx4eZmpqql9vL0mbUpKZpY455SJJDWGgS1JDdB3oSbYneSXJ9zoc25/k/SSvtl+P9bZMSdJKVjOH/ihwAvjFJY6/UFX3XHlJkqTL0dUVepIh4FeBr61tOZKky9XtlMtXgN8ELizT5zNJXkvybJIbO3VIMpZkKsnU/Pz8KkuVpM3t2IFJ5nYMcyHbmNsxzLEDkz09/4qBnuQe4N2qOr5Mt5eBvVV1E/BV4LudOlXVRFWNVNXI4GDHZZSS1EjHDkxy8+Exhs7PsI1i6PwMNx8e62mod3OF/lng15JMA98C7kjyBws7VNUHVXW6vX0UGEiyu2dVStImNzxxkKs4e1HbVZxleOJgz95jxUCvqt+uqqGqGgbuB75fVb++sE+STyVJe/vW9nlP9axKSdrkPn1+dlXtl+Oy16EneSjJQ+3d+4A3krwGPAHcX35yhiT9zNvb96yq/XKsKtCr6vmPliZW1ZGqOtLefrKqbqyqm6rqtqr6Yc8qlKQGmB4b5ww7L2o7w06mx8Z79h4+KSpJ6+D2p0Z55eEJ5rbv5QJhbvteXnl4gtufGu3Ze6RfMyMjIyPlL+eSpNVJcryqRjod8wpdkhrCQJekhjDQJakhDHRJaggDXZIawkCXpIYw0CWpIQx0SWoIA12SGsJAl6SGMNAlqSEMdElqCANdkhrCQJekhjDQJakhug70JNuTvJLkex2OJckTSU4meT3JLb0tU5K0ktVcoT8KnFji2F3ADe3XGHD4CuuSJK1SV4GeZAj4VeBrS3S5F3i6Wl4ErklybY9qlCR1odsr9K8AvwlcWOL4dcBbC/bn2m0XSTKWZCrJ1Pz8/GrqlCStYMVAT3IP8G5VHV+uW4e2Sz6stKomqmqkqkYGBwdXUaYkaSXdXKF/Fvi1JNPAt4A7kvzBoj5zwPUL9oeAt3tSoSSpKysGelX9dlUNVdUwcD/w/ar69UXdngEeaK92uQ14v6re6X25kqSl7Ljcb0zyEEBVHQGOAncDJ4GzwIM9qU6S1LVVBXpVPQ88394+sqC9gC/2sjBJ0ur4pKgkNYSBLkkNYaBLUkMY6JLUEAa6JDWEgS5JDWGgS1JDGOiS1BAGuqTmmpyE4WHYtq31dXKy3xWtqct+9F+SNrTJSc59YYwdH55t7c/MtPYBRkf7Wdma8QpdUiOdfvTgz8O8bceHZzn96ME+VbT2DHRJjbTz1Oyq2pvAQJfUSLPsWVV7Exjokhrp0K5xzrDzorYz7OTQrvE+VbT2DHRJjbTv8VEeGZhgmr1cIEyzl0cGJtj3eDNviIKrXCQ1VGshyyj7D44yOwt79sD4eGMXuAAGuqQGGx1tdoAv5pSLJDXEioGe5BeS/CjJa0neTPLlDn32J3k/yavt12NrU64kaSndTLn8FXBHVZ1OMgAcS/JsVb24qN8LVXVP70uUJHVjxUBvfwD06fbuQPtVa1mUJGn1uppDT7I9yavAu8BzVfVSh26faU/LPJvkxiXOM5ZkKsnU/Pz85VctSbpEV4FeVeer6leAIeDWJL+8qMvLwN6qugn4KvDdJc4zUVUjVTUyODh4+VVLki6xqlUuVfXnwPPA5xa1f1BVp9vbR4GBJLt7VKMkqQvdrHIZTHJNe/vjwJ3Any7q86kkaW/f2j7vqZ5XK0laUjerXK4Ffj/JdlpB/YdV9b0kDwFU1RHgPuDhJOeAvwTub99MlSStk25WubwO3Nyh/ciC7SeBJ3tbmiRpNXxSVJIawkCXpIYw0CWpIQx0SWoIA12SGsJAl6SGMNAlqSEMdElqCANdkhrCQJekhjDQJakhDHRJaggDXZIawkCXpIYw0CWpIQx0SWoIA12SGqKbzxT9hSQ/SvJakjeTfLlDnyR5IsnJJK8nuWVtypUkLaWbzxT9K+COqjqdZAA4luTZqnpxQZ+7gBvar33A4fZXSdI6WfEKvVpOt3cH2q/FHwB9L/B0u++LwDVJru1tqZKk5XQ1h55ke5JXgXeB56rqpUVdrgPeWrA/125bfJ6xJFNJpubn5y+zZElSJ10FelWdr6pfAYaAW5P88qIu6fRtHc4zUVUjVTUyODi46mIlSUtb1SqXqvpz4Hngc4sOzQHXL9gfAt6+ksIkSavTzSqXwSTXtLc/DtwJ/Omibs8AD7RXu9wGvF9V7/S6WEnS0rpZ5XIt8PtJttP6AfCHVfW9JA8BVNUR4ChwN3ASOAs8uEb1SpKWsGKgV9XrwM0d2o8s2C7gi70tTZK0Gj4pKkkNYaBLUkMY6JLUEAa6JDWEgS5JDWGgS1JDGOiS1BAGuiQ1hIEuSQ1hoEtSQxjoktQQBrokNYSBLkkNYaBLUkMY6JLUEAa6JDWEgS5JDWGgS1JDdPMh0dcn+eMkJ5K8meTRDn32J3k/yavt12NrU64kaSndfEj0OeA3qurlJJ8Ajid5rqp+sqjfC1V1T+9LlCR1Y8Ur9Kp6p6pebm//BXACuG6tC5O0AUxOwvAwbNvW+jo52e+KtIxVzaEnGQZuBl7qcPgzSV5L8mySG5f4/rEkU0mm5ufnV1+tpPUzOcm5L4zBzAxUwcxMa99Q37BSVd11TK4G/gQYr6rvLDr2i8CFqjqd5G7g8aq6YbnzjYyM1NTU1GWWLWmtnd49zNWnZi5t37WXq9+bXv+CBECS41U10ulYV1foSQaAbwOTi8McoKo+qKrT7e2jwECS3VdQs6Q+23lqdlXt6r9uVrkE+F3gRFUdWqLPp9r9SHJr+7ynelmopPU1y55Vtav/urlC/yzweeCOBcsS707yUJKH2n3uA95I8hrwBHB/dTuXI2lDOrRrnDPsvKjtDDs5tGu8TxVpJSsuW6yqY0BW6PMk8GSvipLUf/seH+WRB+FL/+8ge5hllj18eWCcOx8f7XdpWkI369AlbUGjowCj7D84yuws7NkD4+MftWsjMtAlLWl01ADfTPxdLpLUEAa6JDWEgS5JDWGgS1JDGOiS1BAGuiQ1hIEuSQ1hoEtSQxjoktQQBrokNYSBLkkNYaBLUkMY6JLUEAa6JDWEgS5JDdHNZ4pen+SPk5xI8maSRzv0SZInkpxM8nqSW9amXEnSUrr5gItzwG9U1ctJPgEcT/JcVf1kQZ+7gBvar33A4fZXSdI6WfEKvareqaqX29t/AZwArlvU7V7g6Wp5EbgmybU9r1aStKRVzaEnGQZuBl5adOg64K0F+3NcGvokGUsylWRqfn5+laVKkpbTdaAnuRr4NvDvquqDxYc7fEtd0lA1UVUjVTUyODi4ukolScvqKtCTDNAK88mq+k6HLnPA9Qv2h4C3r7w8SVK3ulnlEuB3gRNVdWiJbs8AD7RXu9wGvF9V7/SwTknSCrpZ5fJZ4PPAj5O82m77HWAPQFUdAY4CdwMngbPAgz2vVJK0rBUDvaqO0XmOfGGfAr7Yq6IkSavnk6KS1BAGuiQ1xKYK9GMHJpnbMcyFbGNuxzDHDkz2uyRJ2jC6uSm6IRw7MMnNh8e4irMADJ2f4ZOHxzgG3P7UaH+Lk6QNYNNcoQ9PHPxZmH/kKs4yPHGwTxVJ0sayaQL90+dnV9UuSVvNpgn0t7fvWVW71DiTkzA8DNu2tb5Oeg9JF9s0gT49Ns4Zdl7UdoadTI+N96kiaR1NTnLuC2MwMwNVMDPT2jfUtcCmCfTbnxrllYcnmNu+lwuEue17eeXhiSVviLoiRk1y+tGD7Pjw4ntIOz48y+lHvYekn0vrIc/1NzIyUlNTU2ty7sUrYqB1Nb/cDwBpI7uQbWy79BeYcoGwrS70oSL1S5LjVTXS6dimuUJfDVfEqGlm6XyvaKl2bU2NDHRXxKhpDu3qfA/p0C7vIennGhnorohR0+x7fJRHBiaYpnUPaZq9PDIwwb7HnULUzzUy0F0Ro6YZHYU7f2+U/Xun2ZEL7N87zZ2/N8qoea4FNs2j/6tx+1OjHKM1l/7p87O8vX0P02Pj3hDVpjY6igGuZTVylYskNdWWW+UiSVuRgS5JDdHNh0R/Pcm7Sd5Y4vj+JO8nebX9eqz3ZUqSVtLNTdFvAE8CTy/T54WquqcnFUmSLsuKV+hV9QPgp+tQiyTpCvRqDv0zSV5L8mySG5fqlGQsyVSSqfn5+R69tSQJehPoLwN7q+om4KvAd5fqWFUTVTVSVSODg4M9eGtJ0keuONCr6oOqOt3ePgoMJNl9xZVJklbligM9yaeSpL19a/ucp670vJKk1VlxlUuSbwL7gd1J5oAvAQMAVXUEuA94OMk54C+B+6tfj59K0ha2YqBX1b9a4fiTtJY1SpL6yCdFJakhDHRJaggDXZIawkCXpIYw0CWpIQx0SWoIA12SGsJAl6SGMNAlqSEMdElqCANdkhrCQJekhjDQJakhDHRJaggDXZIawkCXpIYw0CWpIVYM9CRfT/JukjeWOJ4kTyQ5meT1JLf0vkxJ0kq6uUL/BvC5ZY7fBdzQfo0Bh6+8LC3n2IFJ5nYMcyHbmNsxzLEDk/0uSdIGsGKgV9UPgJ8u0+Ve4OlqeRG4Jsm1vSpQFzt2YJKbD48xdH6GbRRD52e4+fCYoS6pJ3Po1wFvLdifa7dpDQxPHOQqzl7UdhVnGZ442KeKJG0UvQj0dGirjh2TsSRTSabm5+d78NZbz6fPz66qXdLW0YtAnwOuX7A/BLzdqWNVTVTVSFWNDA4O9uCtt563t+9ZVbukraMXgf4M8EB7tcttwPtV9U4PzqsOpsfGOcPOi9rOsJPpsfE+VSRpo9ixUock3wT2A7uTzAFfAgYAquoIcBS4GzgJnAUeXKtiBbc/NcoxWnPpnz4/y9vb9zA9Ns7tT432uzRJfZaqjtPda25kZKSmpqb68t6StFklOV5VI52O+aSoJDWEga6WyUkYHoZt21pfJzfvunYfvNJWteIcuraAyUnOfWGMHR+217fPzLT2AUY319z8Rw9efbRWf+j8DJ88PMYx8D6DGs85dHF69zBXn5q5tH3XXq5+b3r9C7oCczuGGTp/6Z9lbvtehs5Nr39BUo85h65l7TzV+aGkpdo3Mh+80lZmoItZOj+UtFT7RuaDV9rKDHRxaFfnh5UO7dp8Dyv54JW2MgNd7Ht8lEcGJphmLxcI0+zlkYEJ9j1+5TcR13vFye1PjfLKwxPMbW/9Wea27+WVhye8IaotwZuiAlqrFA8ehNlZ2LMHxsevfIHL4hUn0LpaNmCly7fcTVEDXWvGFSdS77nKRX3hihNpfRnoWjOuOJHWl4GuNeOKE2l9GehaM644kdaXN0UlaRPxpqgkbQEGuiQ1hIEuSQ1hoEtSQxjoktQQfVvlkmQeuPS58P7aDbzX7yI2OMdoeY7Pyhyj5a00PnurarDTgb4F+kaUZGqp5UBqcYyW5/iszDFa3pWMj1MuktQQBrokNYSBfrGJfhewCThGy3N8VuYYLe+yx8c5dElqCK/QJakhDHRJaogtGehJPpfkfyY5meTfdzg+muT19uuHSW7qR539stL4LOj3j5KcT3Lfeta3EXQzRkn2J3k1yZtJ/mS9a+ynLv6N/bUk/yXJa+3xebAfdfZLkq8neTfJG0scT5In2uP3epJbujpxVW2pF7Ad+N/A3wQ+BrwG/P1Fff4x8Mn29l3AS/2ueyONz4J+3weOAvf1u+6NNkbANcBPgD3t/b/e77o32Pj8DvCf2tuDwE+Bj/W79nUco38K3AK8scTxu4FngQC3dZtBW/EK/VbgZFX9WVV9CHwLuHdhh6r6YVX93/bui8DQOtfYTyuOT9u/Bb4NvLuexW0Q3YzRvwa+U1WzAFW1lcapm/Ep4BNJAlxNK9DPrW+Z/VNVP6D1Z17KvcDT1fIicE2Sa1c671YM9OuAtxbsz7XblvJvaP2k3CpWHJ8k1wH/EjiyjnVtJN38HfrbwCeTPJ/keJIH1q26/utmfJ4E/h7wNvBj4NGqurA+5W0Kq80pAHasWTkbVzq0dVy7meSf0Qr029e0oo2lm/H5CvBbVXW+dYG15XQzRjuAfwj8c+DjwH9P8mJV/a+1Lm4D6GZ8/gXwKnAH8LeA55K8UFUfrHFtm0XXObXQVgz0OeD6BftDtK4SLpLkHwBfA+6qqlPrVNtG0M34jADfaof5buDuJOeq6rvrUmH/dTNGc8B7VXUGOJPkB8BNwFYI9G7G50HgP1Zrwvhkkv8D/F3gR+tT4obXVU4tthWnXP4HcEOSX0ryMeB+4JmFHZLsAb4DfH6LXFEttOL4VNUvVdVwVQ0DfwQc2EJhDl2MEfCfgX+SZEeSncA+4MQ619kv3YzPLK3/vZDkbwB/B/izda1yY3sGeKC92uU24P2qemelb9pyV+hVdS7JI8B/o3U3/utV9WaSh9rHjwCPAbuAp9pXoedqi/x2uC7HZ0vrZoyq6kSS/wq8DlwAvlZVHZeoNU2Xf4f+A/CNJD+mNb3wW1W1ZX6lbpJvAvuB3UnmgC8BA/Cz8TlKa6XLSeAsrf/RrHze9hIZSdImtxWnXCSpkQx0SWoIA12SGsJAl6SGMNAlqSEMdElqCANdkhri/wPY2TvDDXs7KAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(xs,ys,c='b')\n",
    "plt.scatter(xs,ys_noise,c='r')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 244,
   "id": "949965ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "SLI=Successive_Linear_Interpolation(xs,ys)\n",
    "SLI_n=Successive_Linear_Interpolation(xs,ys_noise)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 245,
   "id": "42c51b6d",
   "metadata": {},
   "outputs": [],
   "source": [
    "y_n=np.zeros_like(x)\n",
    "y=np.zeros_like(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 246,
   "id": "d9be68a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "for k in range(x.shape[0]):\n",
    "    y[k]=SLI.calculate_I(x[k],order=xs.shape[0]-1,check=False)\n",
    "    y_n[k]=SLI_n.calculate_I(x[k],order=xs.shape[0]-1,check=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 247,
   "id": "f303a2c0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgAAAAFNCAYAAAB2VtfVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABD4UlEQVR4nO3deXhV1b3/8feXDASZwhBmlKCABAIBwqDIYFWwiKJcqVitI6K22N4+V6rW3v5aq61XemuvWge0Fm0tWLUoVVREhYhAmYwQUGSGCEIAmacM6/fH3gkhnCQn5Jyck5zP63n2c4Y9fc/Kydnfvfbaa5lzDhEREYkt9SIdgIiIiNQ8JQAiIiIxSAmAiIhIDFICICIiEoOUAIiIiMQgJQAiIiIxSAmA1Hpm9qyZ/XcE9nu2mR0ys7ia3nddYGadzMyZWfwZrv9zM3sh1HGFav9mdouZLQjDfleb2fAK5s8zswmh3m+4RfrvGYuUAAhmdpGZLTSz/Wa218w+NbP+kY4rWM65u5xzvwnX9sv7QXXObXXONXLOFYZr31VlZsPNLDfIZat1AK5JgT6Xc+63zrmIHehK7/9My9LMEs1st5k1qsJ+ezjn5vnr/8rM/lalwKNUpP+esSjq//ElvMysCfA2cDfwDyARGAIcj2RcUj4zi3fOFSiOOmEokO2cOxTpQKpKf//aTzUA0hXAOTfdOVfonDvqnJvjnFsJp59hlD3TMbPmZvYXM9tuZt+a2Zullh1jZtlmdsDMNpjZ5f77Tc3sz2a2w8y+NrOHi6vRzew8M5vv10bsNrNX/ffNzB43s13+vJVm1tOfN83MHvaff2Fmo0vFEO9vp6//epBf27HPzD6vqCq1MgHKYp6Z/cavQTloZnPMrGWp5cvdt5nd6sd+0Mw2mtmdpeYNN7NcM7vPzL4B/lKFGCuKKct/3OdfyrjAX+c2P5Zvzex9Mzun1Pacmf3IzNYB60q992M/7t1mNsXM6vnz6pnZL8xsi/+3e9nMmpYTa8AyMLOGwLtAOz/OQ2bWLsB38yrzqsf3+Z+7e6l5m83sXv97s9/MXjWzpHLi2GJm/fznN/qfL81/PaH4O15m/wHL0l/u935ZbjKz75bZ3ShgtpldbGarSq0z18yWlHq9wMyuLvVZLjXv/+nnwHX+Pj8vtd1zKvgeVlROzszOK/W69P9Wpd9D8y97lPeZ/b/bLPNqGteb2R2l5pWUp5klmdnfzGyPH+dSM2vtzyv390OqRgmAfAUUmtlLZvZdM2tWxfX/CpwF9ABaAY8DmNkA4GVgMpCMd6az2V/nJaAAOA/oA4wAiqv+fgPMAZoBHYAn/fdH+Nvo6m/vOmBPgHimA9eXej0S2O2cW2Fm7YF3gIeB5sC9wBtmllLFz1yR7wO34pVFor8Pgtj3LmA00MRf/3HzkxZfG3+9c4CJoYgJrzwBkv1LGYv8g8zPgbFACvAJXpmWdjUwEEgr9d41QCbQFxgD3Oa/f4s/XQx0BhoBT5UTZ8AycM4dBr4LbPfjbOSc2156RTPr6sf5n37cs4F/mVliqcW+B1wOpAK9/LgCmQ8M958PBTYCw0q9nh9gndPK0n89EFgLtAQeA/5sZlZqvVF434tFwHlm1tK8hLIn0MHMGptZA6Af3t+ihHPuPeC3wKv+PnuXml3e9zCYcqpIMN/Dij7zdCAXaAdcC/zWzC4JsI2bgaZAR6AFcBdw1J9X0e+HVIESgBjnnDsAXAQ44Hkgz8/QW1e2rpm1xfthvss5961zLt85V/zjeDvwonPuA+dckXPua+fcl/52vwv8p3PusHNuF17SMN5fLx/vx6Wdc+6Yc25BqfcbA+cD5pz7wjm3I0BYfweuMrOz/Nff998DuBGY7Zyb7cf0AbAM70c4VP7inPvKOXcU75JKRjD7ds6945zb4Dzz8ZKgIaW2WwT8P+fccX/boYgpkDuB3/nlW4B3gMmwUrUA/vy9ZeL4H/+9rcAfOZmE3QD8wTm30a/mfgAYbwGulQdRBhW5DnjH/77lA78HGgAXllrmCefcdufcXuBfFZTDfE4e8IcAvyv1ehiBE4DybHHOPe+3E3kJaAsUn8l2BhKcc2udc8fwvg9D8RKplcACYDAwCFjnnAuU8JanvL95MOVUkWC+hwE/s5l1xPutuc//384GXgB+EGAb+XgH/vP8msnlzrkDQfx+SBUoARD8H/tbnHMd8M482uH9iFemI7DXOfdtOfM2BHj/HCAB2OFX7e0DnsM7UwH4GWDAEr+a8jY/xo/wzhz/BOw0s6nmtV8o+1nWA18AV/pJwFWcTADOAcYV79ff90V4P1Ch8k2p50fwzngr3bdf+7LYrxrdh5cYtCy1rTz/IBHKmAI5B/i/UjHuxft7tC+1zLYA65V+bwvedwj/cUuZefH4B8HSgiiDipyyH+dckR9T6biDLYf5wBAzawPEAa8Cg82sE95ZaXaQMZ2yT+fcEf9p8X6vwDsDL73f4ZysZZiHl3BUNek4Zb+c+lmDKaeKBPM9LO8zt8P7vThYatkt5ez7r8D7wAzzLi8+ZmYJVP77IVWgBEBO4Zz7EpiGlwgAHMar4i/WptTzbUBzM0sOsKltwLnlvH8caOmcS/anJs65Hv7+v3HO3eGca4d3Nvp08TVJ59wTzrl+eJcbuuJdXgik+DLAGGCNnxQU7/uvpfab7Jxr6Jx7tNwCCZ1y921m9YE38M7GWjvnkvEODKWrisMxbGegbW4D7iwTZwPn3MJK1utY6vnZQHEV/Xa8H+3S8wqAnaVXDqIMKvv8p+zHr3LuCHxdyXqn8b8vR4AfA1n+AesbvCrvBf5B87TVqrofTlb/FyubABTXRFSUAFR1v5WV0xHK/38/k/2V3XdzM2tc6r2zCfA38msTf+2cS8OrnRgN3EQlvx9SNUoAYpyZnW9m/2VmHfzXHfEOnov9RbKBoebd894UrwoXAL8K/l28g3QzM0sws+JroX8GbjWzS8xrCNbezM7315kD/K+ZNfHnnWtmw/z9jyuOBfgW7wen0Mz6m9lA/yzgMHAMKO/2uxl41wXv5uTZP8Df8GoGRppZnN/QaHip/VUk3l++eEoIYp3SKtp3IlAfyAMKzGs0NaKK2z8TeXhVup1Lvfcs8ICZ9YCSBlfjgtjWZP870BH4Cd5ZM3jJ2E/NLNW8W92Kr1mXbT1eWRnsBFpYOQ0I8aq5r/C/bwnAf+EdKBaWs3xl5gOTOHngnVfmdVmByrJc/nX9Af52iy0EuvnvL3HOrcY7WA/kZCPDsnYCncxvdBmEysopG/i+/x29nJOXPqrNObfN38/v/O9/L7xLha+UXda8RpHp5jXuO4B3SaCwst8PqRolAHIQ7wfm32Z2GO/An4P3w4B/rfpVvGuSy/FuGSztB3j/nF/iNeL6T3+9JfgNuYD9eD+cxWceN+H94K/BO8i/zslq+P5+LIeAWcBPnHOb8BqGPe8vvwWvAeDvA30g/0diEd6Zw6ul3t+GVyvwc7wf7G14tQjB/B88g9cIqXgKuiV+Zfv2zzB/jPfj/C1eu4VZVdn+mfCrZx8BPvWrUwc552YC/4NX9XoA77tQtuV6IG/hfT+y8c5q/+y//yJedW4WsAkvcbsnQCwVloFfMzUd2OjH2q7M+mvx2lk8CewGrgSudM6dCCL2QObjtTnJKud12fhPK8tKtn8JsKh0dbrzGjuuAFaXinsR3jX1XeVs5zX/cY+Zrahkn8GU00/89/bhtd94s7JtVtH1QCe82oCZeO0JPgiwXBu834UDeJf05uMl0VDx74dUgTkXjppFEYkVZuaALqUutUglzOxpIMc593SkY5HYpY6ARERqXjbenQgiEaNLACKAnexgpuwU7G1oNca8PtMDxfpupGOT4DjnprrAt7GK1BhdAhAREYlBqgEQERGJQUoAREREYlBMNQJs2bKl69SpU6TDEBERqRHLly/f7ZwLON5JTCUAnTp1YtmyZZEOQ0REpEaY2Zby5ukSgIiISAxSAiAiIhKDojYBMLPNZrbKzLLN7LR6e/M8YWbrzWylnTp2uoiIiFQg2tsAXOyc213OvO8CXfxpIF5f7QNrKjAREZHaLGprAIIwBnjZeRYDyWamASFERESCEM0JgAPmmNlyM5sYYH57vBHViuX6753CzCaa2TIzW5aXlxemUEVERGqXaE4ABjvn+uJV9f+o1DjzxSzAOqf1a+z3uZ3pnMtMSQl4K6SIiEjMidoEwDm33X/chTdu9IAyi+QCHUu97oA3xrSIiIhUIioTADNraGaNi58DI4CcMovNAm7y7wYYBOzX6FoiIiLBida7AFoDM80MvBj/7px7z8zuAnDOPQvMBkYB64EjwK0RilVERCQEVgBbgatrZG9RmQA45zYCvQO8/2yp5w74UU3GJSIiEnoO+BPwX0BnYDQ1cXiOyksAIiIisSA/Pw+4FrgHDg6CxVOoqXNzJQAiIiI1zO3ezYqnbuTQtla4/H/iflEfmmbBhPtrLIaovAQgIiJSZzgHX30FCxbAwoWwbR5260b6TgK3Bb4aCl0zboW/DIJBg2osLCUAIiIioVRQAJ99BllZ3kF/wQLYvRtaAQ8nwXPHoSiB5f/swhW3rOGmuybz2GOP1XiYSgBERESqIz8fli6FefO8g/6nn8KhQ9688zrD5D7wH/sh9TOwfLC7gP+m39i2fDM2cmErARAREamKggJYsQI++gg+/tg7wz9yxJvXsyfc/AO4qiMM2gNN3gA+AJoDdwOT8MawizwlACIiIhVxDtasgblzvYP+vHlw4IA3r2dPuO02GNULBgNNFgL/BHbi9Vh/MfA7vHv7kyIRfbmUAIiIiJT1zTfwwQfeNHcu7PA7mj33XLjlahjTCfrHQePVwJvAU/6KzYCReP3UjcS78B+dlACIiIgcP+5V5b//PsyZA59/7h0hM5vC5DQY1g+6nYCGXwEvl1rxbOAi4AJ/6kNtObTWjihFRERCbdMmePddeO9dWPchnHcUeteD37eE9DaQsgfq7QcW4XWb0wUYCNwJZPhTm0hFX21KAEREJDacOAEL5sHyv8HeD6DNN94x/IZ60LTIX6gI71p9OtDDf+wJnE+0XcOvLiUAIiJSJ+XnnyBh3ypY/QIc/gBaboJBRfAdf4ETiVCYBkkD8TKBXngH+yaRCrlGKQEQEZE6oghcDu7r6ezf9jfOOjsX2gPDgaMGua1hR39oPwaSBkNiFyAusiFHkBIAERGppRywHgrehT3/gEbLoOFxrINXo79nPixaCKPveQXrei10SYx0wFFFCYCIiNQiB4C5cOItyJ8NDXd7R7KjwPtx8G0faDue38zfzP97+hkmT57Mld2+H+GYo5M55yIdQ43JzMx0y5Yti3QYIiJSJWuBWXD8TYhfDHFFsB+YC/y7EcRdDhfcCJdeBmedFdlQo4yZLXfOZQaapxoAERGJMkV4t969BSfegMSN3ttfAu8Cn7eDjt+DMdfC7wZBXOxex68OJQAiIhIF8oH54N6AgtchYTfkG8xz8BawMQ0uHA9XXw339QSzCMdb+ykBEBGRCMkHPgL3Dyh8A+L3e63133EwE9g7EEaMg3uvgdTUSAdb5ygBEBGRGlQAzAM3Awpeg4QDcMjgLQf/rAfHh8Lo78EfroE2tbeXvdpACYCIiIRZEbAQ3HTInw6J35Y66MdB0aVw1ffg+THQokWkg40ZSgBERCQMHPA5uFfg+MuQtAuOAf8C3oiHwpEw5jp48UpITo5sqDFKCYCIiITQBij6Oxz9MzTc4tX4zwXeSIDjI2H09fD8aGgSG93tRjMlACIiUk27oGg6HJoKTdZ4A+ctA16Lh4Mj4bs3whNXQOPGkQ5USlECICIiZ+AgFP0T9j0NTZdCnIONeAf9vEvhOzfDo6OhUaNIByrlUAIgIiJBOgFFs2HPk9A0CxILvB75no+DHRfDBRPggSt00K8lojIBMLOOwMtAG7zmo1Odc/9XZpnheN1DbPLf+qdz7qEaDFNEJAYUQdF82PU4NJ4DDY97b78YB1uHQa874UdX6qBfC0VlAoDXbOS/nHMrzKwxsNzMPnDOrSmz3CfOudERiE9EpA5zULQMdvwBGr4NyYegETCrHmy4AM67G264Wtf0a7moTACcczuAHf7zg2b2Bd6ozmUTABERCYH8/HwS4tbB17+H+m9Cq28hBZhTD9ZmQoe7YfQ4HfTrkKhMAEozs05AH+DfAWZfYGafA9uBe51zq2syNhGR2s4Vrmf9R3eRnPIhKRl4p1pZBq/1gtY/ghHjYbRu2auLojoBMLNGwBvAfzrnDpSZvQI4xzl3yMxGAW8CXQJsYyIwEeDss88Ob8AiIrVB4SbY/HuIfx07ZxddLoOiRbDyx5A+7EnssptguA76dV29SAdQHjNLwDv4v+Kc+2fZ+c65A865Q/7z2UCCmbUMsNxU51ymcy4zJSUl7HGLiESlgo2w7m7Y2AbiOsO5T8PuPHipB2/ddjHJF8LfkiZj/zFJnfTEiKisATAzA/4MfOGc+0M5y7QBdjrnnJkNwEtm9tRgmCIi0e3EWtg4BRJnQec8r4402+CVnpA8AYbdDv0aMeZmOPBipIOVmhaVCQAwGPgBsMrMsv33fg6cDeCcexa4FrjbzAqAo8B455yLQKwiItHjWDZs/D00eBdS98L5QHY9mNELmk2EIbdCxlmRjlKiQFQmAM65BYBVssxTwFM1E5GISLRycHAebPkjJH8MHQ5CGrA0DpZkQsqdcOGNkJEU2TAl6kRlAlArOAf5+XD0KBw7dvIxP9+bCgq8qajIm5zzJjOoV+/klJBwckpMhKQkb2rQwHuMi4v0JxWRqHMCdr8O3zwHrf8NKce9M/1FifDpEGh3Nwy6FvonRDpQiWJKAM7UL34Bv/1t+PdTvz40bOj1stWokXcPbtOm3tSkCTRrBs2bn3xs0QJatoSUFO95/frhj1FEws/tga9fgAOvwNlroGUhNAAWNoS934FzfwKDL/NOLESCoATgTF12mXdgLn22Xr++dxafkADx8d4UF3fyrN/MqwUorhUoLDxZY5CfDydOwPHjXm1C8XTkCBw6BIcPw8GD3rR/P2zd6j3u2+ctV54mTaB165NTmzbQrp03tW3rPXbo4I3HbYGvuuTn55OQoDMJkZrlIH8VbH0G3NvQKRc6AN8Ac1vCicsh/adwaZ9y/3dFKmKx1G4uMzPTLVu2LNJhhN6xY/Dtt7B3L+ze7U15ed7jrl2wc6c37doFO3Z4SUNZZzWAjLbQPxk6tIA2LXBtWpCzazu//8tMbvjJr7nsipsxSwHUgEgkPI7Cwbdh+wvQ7FNoddh7e6XBms5Q/1oY8CNo3zGyYUqtYWbLnXOZAecpAYhBR4/Cji1w/E1I/AAabIDkb+Cs48Gtf6A+7E+BE+dAfHdoMhCSh4GdC9RTjYFIVbh1FOY+T9yJt6D9ekgqgiPAgkT4ph+0ugUGX68ueOWMVJQA6BJATHFAFjR4BTq/DnwLNAN6AaOA7kBXoBHkA3n7+Mf/PcVXn77HFRd1o0/X5lCwAxrlQZtc6JILHT8FXvC2frgeJzacRU7WIZIbDaNzz1uxDpdC63aqohQpcRCOvA3fvIRrtBBrdZC4jsA6cDNTsPxLoNudcOlQXc+XsFINQMxYDfwQyMKrwr8GuAG4FDiDs/Vjx2DLFti0Gg4shsIVuIYbONp0Mw36gBV3JHYYyI6DTS1hbxdw/aHtQDi/O3Tp4rWfEKnTCqDo37DjZSh6H9puhXgHh8BlxfP57AL+8R5M2wBfFxVhSpYlhFQDENMOA78B/hdoAjyD18dSw+ptNikJunXzJsYCXscNv/rZz3j84in86cFrmTiuExQtgXZroX8eJO4EFsBOYAnwGrClNRztCR3SoXt3OP9871HdNkut5YDVsPc1ODgLUlbDWfnQFlgOzGkFhZdA99uxy4by93kPMmXDFCZPnqyDv9Qo1QDUaZ/gHey3ALcC/4M3vmcknABWwolP4NDHEL8cmmw/OXuDwVIHy/Cmrc2gQ08vGUhLO/nYvr0uJ0iUccBaOPAv2DcTkrOhiX9nznpgUQPY3x9aXw/DxkKrVhGMVWKNGgH6YisBeA+vmr8j3rAKQyIbTkD78Y72S8Atg8J/Q/zXJ2d/nQTLC+Hf+fA5sBLY1wi6p0GPHl5CUDydfbaul0oNKQLWwMF3YN8saJINTY94s3KBT+Phmx7Q5GoY+D0veVXSKhGiBMAXOwnATOA6oCfwPpE76z8TeXj1pJ8BK8B9Brbh5Owj9WF9Eqw4AcuOek0bVgOHG3iJQVqZKTVVvSlKNZ3wEtRv/wWH50CzNdDomDcrF1gYB9u7QYNR0Oda6NvP6wNEJAooAfDFRgLwCnAzMACYDSRHNJrQOADk4FUBFE+rgX0nFzmcBJuTIOcELD8CXwHrgNxESO1++qWE887zOm0SOc1OyM+CvLfALYSULZBY5M36CliSAHnd4azLoffV0C/T6/xLJAopAfDV/QTgr3gH/+HALKBRRKMJLwfs4GQVwFrgS3/65uRiRcDuJO9a7OpjsBHYAGytB6RCmx5ezUFxA8Ru3bxuliVGHIWi5bB7Nhz5GBqvgRYHvFnHgRVATmM4lAHJo6DfKO/yk2qVpJZQAuCr2wnAGiATGAS8g9dJeKzah3f6v95/XAdsALceLO/URQ/Xg01FsAmvreQW4EAyxHWGxj2gbQZ06+4lBuecox/+Wu0IFGVD3vtweD40WA2t9kCc/xu4BVgeBzs6AhdC+ythwFCvu2yRWkoJgK/uJgDHgYHAdrzq8TaRDSeqHcKrBtjkT5uhaAOcWAdx2yDhyKmLnwC2AVuBrw0ON4eidpB4LjTpAa36Qqe+5LduQ0KShluNDn7t0IEFsPtjKFgGjdZD631QnL/lASsMctvAiV7Q9DJIH+HVBOn6vdQh6gegznsAr5n82+jgX5lGeD0f9jr5Vj2g5Ni9n5KqALcVCtZC4y+g+2bI2AVN9kLcHmAV8Ka3ylEvdyjYlUjc8ZZYYRuIPwfO6gLN06F1X2jSFf27hZoDtwN2fwJ7F8CJzyFpI7TeCU0KvG4vmuAlb9kJkHcOFKRD42Fw/mXwnTRdu5eYpl+kWu894HFgEnBFhGOpC5pSkiAYXqeJp4x9VAh8A0VbYO8q3N6VzHnrac7uCCntT9Ci625I2Q71V5y62UJgbzzsbwhHm0Nha4jvAA06Q5Nu0CwN4jsCrSivZ8ZYGGPh9M/o4PhW2LUY9q+A42ug3gZovAPaHIBGRd5NLil4V37W1oNVLeBIZ4jvCy2Hw/kXwnfVf4RIWUoAarVdwC14t/s9FtlQYkYc0B7qtYeWF2It4cO8hkz5mdeT22OPPQauEPaug2+Ww7c5cPgrKNgKcTvhrG+h6RZovynw3ZlFwIFEONQIjjeFwhY4a8P6LUd58e8fMnr8JC689HqsXmugJd4pbi09sLkiOLAV9qzCHfiC9aveZ92XH3HhRR1p2uYElrwfWh+Hs5zXnUVHvPLZCuQ2hC/O9gakSugNzS+CzhfCAI07IRIstQGo1cYB/8LrTKdnhGORKtm3D77eBLtWwf4v4cgGKPwabBfU3wuNDkLyCa9CoDXl39BRYHAowesH4XhDyG8ERY3BNQVrBpYMccmQ0Nyb4ptCQjIk+lNCI7AkoD5Qn/x8R0JCfSpMKoryIf8onDgMx/bD8X1w7FvIPwDH9sCJPZC/Bwq+haK94PZC3H6IPwhJh6DhUUg+Di2KIMCdmEW7wHYmYfsaw5FWUHQOJJ4PTfpA6wHQ/lw1xhQJktoA1EnLgdeBX6GDfy2UnAzJfaBHn/KXKSiAvDxYv5Npv/81Kz99k2su78uQgeeC2wX19kD8Aah/CBoc9ZKGxkVe1w/FUxU7RyyufXcFYIV47ems1BTnb7O+PwUzQm0RsN/gQDwcToSDTWB3E1jbHGjlt5c4j+dfXcB/P/E6N026l8emTKla4CJSZaoBqLWuABbjtWRvUsmyEjMKCuDgQdi/Hw7sh6O74fhu74w8fy8UHYSiQ+AOgzsEdgI4jrMTfLpwHvXioV4cDByUidUr8qrTzfCO+gb14oEEsESwBKjXAOo19Kb4xl4NQ1KKN53VBhq2hYbtwXTGLhIJqgGocxbh9fL3O3Twl1PEx0OzZt5UBQbMyvoZU37ntWUYNFJtSkTqOtUA1EqX4d32t5G63dufiIhUh2oA6pT5wFzgD+jgLyIiZ0rjp9YqDvhvoB1wV4RjERGR2kw1ALXKXOAT4Cliu69/ERGprqitATCzy81srZmtN7P7A8w3M3vCn7/SzPpGIs6a9RBebygTIh2IiIjUclGZAJhZHPAn4LtAGnC9maWVWey7QBd/mgg8U6NB1rgvgQXAj/FuwBYRETlzUZkAAAOA9c65jc65E8AMYEyZZcYALzvPYiDZzNrWdKA15yW8XlhujHQgIiJSB0RrAtAebxDWYrn+e1Vdpo4oBF7Gq/TQaH8iIlJ90ZoABOqIvGyHBcEsg5lNNLNlZrYsLy8vJMHVvDnAduDWSAciIiJ1RLQmALl4rd2KdcA7AlZ1GZxzU51zmc65zJSUQMOv1QbTgBbA6AjHISIidUW0JgBLgS5mlmpmicB4YFaZZWYBN/l3AwwC9jvndtR0oOG3F3gTuIGAQ6eJiIicgajsB8A5V2Bmk4D38Vq+veicW21md/nzn8XrDH8UsB44Qp2tH58BnKDOfjwREYkIjQUQ9foD+UB2hOMQEZHapqKxAKL1EoAAkAMsQ2f/IiISakoAoto0IAHv+r+IiEjoKAGIWkXAK8AVQMsIxyIiInWNEoColQ18A4yNcBwiIlIXKQGIWu/6jyMjGoWIiNRNSgCi1rtAJtAq0oGIiEgdpAQgKn0LLAIuj3QgIiJSRykBiEof4DUC/G6kAxERkTpKCUBUeg9oBgyMdCAiIlJHKQGIOg4vARiB1wuyiIhI6CkBiDqfAzvQ9X8REQknJQBRp/j2PyUAIiISPkoAos67QB+gTaQDERGROkwJQFTZDyxEZ/8iIhJuSgCiylygEN3+JyIi4aYEIKq8CzQFLoh0ICIiUscpAYgaxbf/XQbERzgWERGp65QARI1NwNfAJZEOREREYoASgKix2H9U9b+IiISfEoCosRg4C+gR6UBERCQGKAGIGouB/uj6v4iI1AQlAFHhGJANDIpwHCIiEiuUAESFz4B8lACIiEhNUQIQFYobAGr4XxERqRlKAKLCYuAcoG2kAxERkRihBCAqLEZn/yIiUpOiLgEwsylm9qWZrTSzmWaWXM5ym81slZllm9myGg4zhHYAW9H1fxERqUlRlwAAHwA9nXO9gK+ABypY9mLnXIZzLrNmQguHf/uPSgBERKTmRF0C4Jyb45wr8F8uBjpEMp7wWwwkAH0iHYiIiMSQqEsAyrgNb4i8QBwwx8yWm9nEGowpxBbjHfyTIh2IiIjEkIh0O2dmc4E2AWY96Jx7y1/mQaAAeKWczQx2zm03s1bAB2b2pXMuK8C+JgITAc4+++yQxB86BcBS4PZIByIiIjEmIgmAc+7Siuab2c3AaOAS55wrZxvb/cddZjYTGACclgA456YCUwEyMzMDbityVgNH0PV/ERGpaVF3CcDMLgfuA65yzh0pZ5mGZta4+DkwAsipuShDpbgDICUAIiJSs6IuAQCeAhrjVetnm9mzAGbWzsxm+8u0BhaY2efAEuAd59x7kQm3OhYDKUBqpAMREZEYE3VDzznnzivn/e3AKP/5RqB3TcYVHovxzv4t0oGIiEiMicYagBixD/gS9QAoIiKRoAQgYlb5j30jGoWIiMQmJQARs9p/7BHRKEREJDYpAYiYHLy2jh0jHYiIiMQgJQARsxrv7F8NAEVEpOYpAYiYHKBnpIMQEZEYpQQgInYBu9H1fxERiRQlABFR3GmhagBERCQylABEhO4AEBGRyFICEBE5QDMCD4goIiISfkoAImI1XvW/7gAQEZHIUAJQ4xxeDYCq/0VEJHKUANS47cB+1ABQREQiSQlAjSu+A0A1ACIiEjlKAGqc7gAQEZHIUwJQ43KAVkBKpAMREZEYVmkCYGaTzKxZTQQTG4rvABAREYmcYGoA2gBLzewfZna5menetTNWBKxB1f8iIhJplSYAzrlfAF2APwO3AOvM7Ldmdm6YY6uDtgKHUA2AiIhEWlBtAJxzDvjGnwrwurF73cweC2NsdZAaAIqISHSIr2wBM/sxcDPe8HUvAJOdc/lmVg9YB/wsvCHWJboFUEREokOlCQDQEhjrnNtS+k3nXJGZjQ5PWHXVaqA9kBzhOEREJNZVmgA4535ZwbwvQhtOXZeDrv+LiEg0UD8ANaYQ+AJV/4uISDRQAlBjNgLHUA2AiIhEAyUANWad/9gtolGIiIhAFCYAZvYrM/vazLL9aVQ5y11uZmvNbL2Z3V/TcVbdBv9R3SeIiEjkBXMXQCQ87pz7fXkzzSwO+BNwGZCL11PhLOfcmpoKsOo2AmfhjQMgIiISWVFXAxCkAcB659xG59wJYAYwJsIxVWIj0BlQT8oiIhJ50ZoATDKzlWb2YjkDEbUHtpV6neu/F8U2oOp/ERGJFhFJAMxsrpnlBJjGAM/gHSkzgB3A/wbaRID3XDn7mmhmy8xsWV5eXqg+QhU5TtYAiIiIRF5E2gA45y4NZjkzex54O8CsXKBjqdcdgO3l7GsqMBUgMzMzYJIQfjuBoygBEBGRaBF1lwDMrG2pl9dwsgP90pYCXcws1cwSgfHArJqI78zoDgAREYku0XgXwGNmloFXb74ZuBPAzNoBLzjnRjnnCsxsEvA+EAe86JxbXc72osBG/1E1ACIiEh2iLgFwzv2gnPe3A6NKvZ4NzK6puKpnI16zhU4RjkNERMQTdZcA6qYNeM0U6kc6EBEREUAJQA3RHQAiIhJdlADUCCUAIiISXZQAhN0RvO4MdAeAiIhEDyUAYbfJf1QNgIiIRA8lAGGnWwBFRCT6KAEIu+IEQJcAREQkeigBCLsNQGOgRaQDERERKaEEIOw0DLCIiEQfJQBhtxFV/4uISLRRAhBWRagPABERiUZKAMJqB3AcJQAiIhJtlACEle4AEBGR6KQEIKw2+I+qARARkeiiBCCsNuIV8dmRDkREROQUSgDCaiPewT8x0oGIiIicQglAWOkOABERiU5KAMJqA0oAREQkGikBCJtDwC50B4CIiEQjJQBho2GARUQkeikBCJviPgBSIxqFiIhIIEoAwibXf9QtgCIiEn2UAITN10ACkBLpQERERE6jBCBscoG2qIhFRCQa6egUNl8DHSIdhIiISEBKAMImF2gf6SBEREQCUgIQFg6vBkAJgIiIRKf4SAdQlpm9CnTzXyYD+5xzGQGW2wwcBAqBAudcZg2FGIQDwGF0CUBERKJV1CUAzrnrip+b2f8C+ytY/GLn3O7wR1VVxbcAqgZARESiU9QlAMXMzIDvAd+JdCxV97X/qBoAERGJTtHcBmAIsNM5t66c+Q6YY2bLzWxieRsxs4lmtszMluXl5YUl0NOpBkBERKJbRGoAzGwu0CbArAedc2/5z68HplewmcHOue1m1gr4wMy+dM5llV3IOTcVmAqQmZnpqhl6kIprANrVzO5ERESqKCIJgHPu0ormm1k8MBboV8E2tvuPu8xsJjAAOC0BiIyv8XoArB/pQERERAKK1ksAlwJfOudyA800s4Zm1rj4OTACyKnB+CqhPgBERCS6RWsCMJ4y1f9m1s7MZvsvWwMLzOxzYAnwjnPuvRqOsQLqA0BERKJbVN4F4Jy7JcB724FR/vONQO8aDqsKvgYGRjoIERGRckVrDUAtdhzIQzUAIiISzZQAhNx2/1F9AIiISPRSAhBy6gNARESinxKAkCvuA0AJgIiIRC8lACGnboBFRCT6KQEIuVygIdAk0oGIiIiUSwlAyH2Nd/ZvkQ5ERESkXEoAQk6dAImISPRTAhBy6gZYRESinxKAkCrC6wdADQBFRCS6KQEIqV1AAaoBEBGRaKcEIKTUB4CIiNQOSgBCSn0AiIhI7aAEIKTUDbCIiNQOSgBC6mu8EZZbRToQERGRCikBCKlcoC0QF+lAREREKqQEIKTUCZCIiNQOSgBCqrgbYBERkeimBCCk1AugiIjUDkoAQuYAcAjVAIiISG2gBCBk1AmQiIjUHkoAQkZ9AIiISO2hBCBkVAMgIiK1hxKAkPnGf2wb0ShERESCoQQgZHYBjYCzIh2IiIhIpZQAhMwu1AWwiIjUFhFJAMxsnJmtNrMiM8ssM+8BM1tvZmvNbGQ56zc3sw/MbJ3/2KxmIq+IEgAREak9IlUDkAOMBbJKv2lmacB4oAdwOfC0mQXqWP9+4EPnXBfgQ/91hCkBEBGR2iM+Ejt1zn0BYGZlZ40BZjjnjgObzGw9MABYFGC54f7zl4B5wH1hCjdIu4D+kQ1BRGq1/Px8cnNzOXbsWKRDkVomKSmJDh06kJCQEPQ6EUkAKtAeWFzqdXl967Z2zu0AcM7tMLMIn3oXAXmoBkBEqiM3N5fGjRvTqVOnQCdIIgE559izZw+5ubmkpqYGvV7YLgGY2VwzywkwjalotQDvuWrGMdHMlpnZsry8vOpsqgL7gAKUAIhIdRw7dowWLVro4C9VYma0aNGiyjVHYasBcM5degar5QIdS73uAGwPsNxOM2vrn/23xat/Ly+OqcBUgMzMzGolE+Ur3r0SABGpHh385Uycyfcm2m4DnAWMN7P6ZpYKdAGWlLPczf7zm4G3aii+cigBEJG64ZFHHqFHjx706tWLjIwM/v3vf4d1f8OHD2fZsmXV2sa0adNISUkhIyODtLQ0nn/++QqXnzBhAmvWrKlwmTfffLPSZWq7SN0GeI2Z5QIXAO+Y2fsAzrnVwD+ANcB7wI+cc4X+Oi+UumXwUeAyM1sHXOa/jiAlACISem3atMHMQja1adOmwv0tWrSIt99+mxUrVrBy5Urmzp1Lx44dK1wnWlx33XVkZ2czb948fv7zn7Nz585yl33hhRdIS0urcHtKAMLEOTfTOdfBOVffOdfaOTey1LxHnHPnOue6OefeLfX+BOfcMv/5HufcJc65Lv7j3kh8jpOUAIhI6FV0EAvH9nbs2EHLli2pX78+AC1btqRdu3YAPPTQQ/Tv35+ePXsyceJEnPOuqA4fPpyf/vSnDB06lO7du7N06VLGjh1Lly5d+MUvfgHA5s2bOf/887n55pvp1asX1157LUeOHDlt/3PmzOGCCy6gb9++jBs3jkOHDgFw//33k5aWRq9evbj33nsr/AytWrXi3HPPZcuWLXz44Yf06dOH9PR0brvtNo4fP14Sc3GtQ6NGjXjwwQfp3bs3gwYNYufOnSxcuJBZs2YxefJkMjIy2LBhA0888URJDOPHjw+2yKNatF0CqKV24bVfbBHpQEREztiIESPYtm0bXbt25Yc//CHz588vmTdp0iSWLl1KTk4OR48e5e233y6Zl5iYSFZWFnfddRdjxozhT3/6Ezk5OUybNo09e/YAsHbtWiZOnMjKlStp0qQJTz/99Cn73r17Nw8//DBz585lxYoVZGZm8oc//IG9e/cyc+ZMVq9ezcqVK0uSivJs3LiRjRs30qFDB2655RZeffVVVq1aRUFBAc8888xpyx8+fJhBgwbx+eefM3ToUJ5//nkuvPBCrrrqKqZMmUJ2djbnnnsujz76KJ999hkrV67k2WefrU4xRw0lACGxE+/gH213VYqIBK9Ro0YsX76cqVOnkpKSwnXXXce0adMA+Pjjjxk4cCDp6el89NFHrF69umS9q666CoD09HR69OhB27ZtqV+/Pp07d2bbtm0AdOzYkcGDBwNw4403smDBglP2vXjxYtasWcPgwYPJyMjgpZdeYsuWLTRp0oSkpCQmTJjAP//5T846K/B4K6+++ioZGRlcf/31PPfcc+Tl5ZGamkrXrl0BuPnmm8nKyjptvcTEREaPHg1Av3792Lx5c8Dt9+rVixtuuIG//e1vxMfXjd/6uvEpIk69AIpI3RAXF8fw4cMZPnw46enpvPTSS4wfP54f/vCHLFu2jI4dO/KrX/3qlFvOii8Z1KtXr+R58euCggLg9FbqZV8757jsssuYPn36aTEtWbKEDz/8kBkzZvDUU0/x0UcfnbbMddddx1NPPVXyOjs7O6jPm5CQUBJLXFxcSbxlvfPOO2RlZTFr1ix+85vfsHr16lqfCKgGICSUAIhI7bd27VrWrVtX8jo7O5tzzjmn5GDfsmVLDh06xOuvv17lbW/dupVFi7xOXadPn85FF110yvxBgwbx6aefsn79egCOHDnCV199xaFDh9i/fz+jRo3ij3/8Y9AH9vPPP5/NmzeXbO+vf/0rw4YNCzrexo0bc/DgQQCKiorYtm0bF198MY899hj79u0raZ9Qm9Xu9CVq7AL6RDoIEZFqOXToEPfccw/79u0jPj6e8847j6lTp5KcnMwdd9xBeno6nTp1on//qnd73r17d1566SXuvPNOunTpwt13333K/JSUFKZNm8b1119f0ljv4YcfpnHjxowZM4Zjx47hnOPxxx8Pan9JSUn85S9/Ydy4cRQUFNC/f3/uuuuuoOMdP348d9xxB0888QQzZszg9ttvZ//+/Tjn+OlPf0pycnLQ24pWVtySMxZkZma66t5vGlgz4AfAE2HYtojEii+++ILu3buXvG7Tpk1I7wRo3bo133zzTci2F6zNmzczevRocnJyanzfsaTs9wfAzJY75zIDLa8agGo7gdcVsC4BiEhoReJgLbFDbQCqrXh8ASUAIiKBdOrUSWf/UUgJQLWpEyAREal9lABUmxIAERGpfZQAVJsSABERqX2UAFSbEgAREal9lABU2y6gPtA40oGIiIgETQlAtRX3AmiVLSgiUjVt2oBZ6KZKhgMGrzvcjIwMevbsybhx4wKO2hesW265paTXwAkTJlQ4vO68efNYuHDhGe8LYNasWTz6qDc6fNnhfEuPAFhdv/3tb0OyndJ++ctfMnfu3JBvtyJKAKpN3QCLSJiEeDjgYLbXoEEDsrOzycnJITEx8bSR7woLC89o1y+88AJpaWnlzg9FAnDVVVdx//33A6cnAKEUjgTgoYce4tJLLw35diuiBKDalACISN00ZMgQ1q9fz7x587j44ov5/ve/T3p6OoWFhUyePJn+/fvTq1cvnnvuOcAb0GfSpEmkpaVxxRVXsGvXrpJtlT4Df++99+jbty+9e/fmkksuYfPmzTz77LM8/vjjZGRk8Mknn5wWS2FhIZ07d8Y5x759+6hXr17J6H7FcU6bNo1JkyaxcOFCZs2axeTJk8nIyGDDhg0AvPbaawwYMICuXbuW7OPYsWPceuutpKen06dPHz7++GOAkm0VGz16NPPmzeP+++/n6NGjZGRkcMMNNwQst82bN9O9e3fuuOMOevTowYgRIzh69Cjgja8waNAgevXqxTXXXMO3334LnFpbcv/995OWlkavXr249957AcjLy+M//uM/6N+/P/379+fTTz89kz/pKdQTYLXtAnpEOggRkZAqKCjg3Xff5fLLLwe8EflycnJITU1l6tSpNG3alKVLl3L8+HEGDx7MiBEj+Oyzz1i7di2rVq1i586dpKWlcdttt52y3by8PO644w6ysrJITU1l7969NG/enLvuuotGjRqVHPDKiouLo2vXrqxZs4ZNmzbRr18/PvnkEwYOHEhubi7nnXdeyRDDF154IVdddRWjR4/m2muvPeUzLVmyhNmzZ/PrX/+auXPn8qc//QmAVatW8eWXXzJixAi++uqrcsvl0Ucf5amnnqp0UKJ169Yxffp0nn/+eb73ve/xxhtvcOONN3LTTTfx5JNPMmzYMH75y1/y61//mj/+8Y8l6+3du5eZM2fy5ZdfYmbs27cPgJ/85Cf89Kc/5aKLLmLr1q2MHDmSL774osIYKqMEoFocsBPVAIhIXVF8dgvemfXtt9/OwoULGTBgAKmpqQDMmTOHlStXlpyx7t+/n3Xr1pGVlcX1119PXFwc7dq14zvf+c5p21+8eDFDhw4t2Vbz5s2Djm3IkCFkZWWxadMmHnjgAZ5//nmGDRsW9OBEY8eOBaBfv35s3rwZgAULFnDPPfcA3giC55xzToUJQLBSU1NLyrF4f/v372ffvn0loxLefPPNjBs37pT1mjRpQlJSEhMmTOCKK65g9OjRAMydO/eUSxoHDhzg4MGDNG585g3QlQBUy0HgONA60oGIiIREcRuAsho2bFjy3DnHk08+yciRI09ZZvbs2ZhV3CDaOVfpMuUZMmQIzz77LNu3b+ehhx5iypQpzJs3j6FDhwa1fv369QGvNqGgoKAknkDi4+MpKioqeV08JHKwivdVvL/iSwCViY+PZ8mSJXz44YfMmDGDp556io8++oiioiIWLVpEgwYNqhRHRdQGoFrUB4CIxJ6RI0fyzDPPkJ+fD8BXX33F4cOHGTp0KDNmzKCwsJAdO3aUXE8v7YILLmD+/Pls2rQJ8Kq8ARo3bszBgwcr3O/AgQNZuHAh9erVIykpiYyMDJ577jmGDBly2rLBbA9g6NChvPLKKyWfY+vWrXTr1o1OnTqRnZ1NUVER27ZtY8mSJSXrJCQklHz2qmjatCnNmjUraX/w17/+taQ2oNihQ4fYv38/o0aN4o9//GNJMjZixAieeuqpkuUquwQRDCUA1aIEQETCqHWIaxdDtL0JEyaQlpZG37596dmzJ3feeScFBQVcc801dOnShfT0dO6+++7TDm4AKSkpTJ06lbFjx9K7d2+uu+46AK688kpmzpxZbiNA8M6qO3bsyKBBgwCvRuDgwYOkp6eftuz48eOZMmUKffr0KWkEGMgPf/hDCgsLSU9P57rrrmPatGnUr1+fwYMHk5qaSnp6Ovfeey99+/YtWWfixIn06tWr3EaAFXnppZeYPHkyvXr1Ijs7m1/+8penzD948CCjR4+mV69eDBs2jMcffxyAJ554gmXLltGrVy/S0tJOuzvjTFh51R91UWZmpgvVfaCeN4FrgBVAnxBuV0RiUaDx3EWCFej7Y2bLnXOZgZZXDUC1qAZARERqJzUCrJbiBCAlolGIiNQljzzyCK+99top740bN44HH3wwQhEFtmfPHi655JLT3v/www9p0aJFBCKqGiUA1bILSAYSIxyHiEjd8eCDD0bdwT6QFi1ahKQxXqToEkC1qBdAERGpnSKSAJjZODNbbWZFZpZZ6v3LzGy5ma3yH0/vRcJb7ldm9rWZZfvTqJqLvjQlACIiUjtF6hJADjAWeK7M+7uBK51z282sJ/A+0L6cbTzunPt9GGMMwi6gW2RDEBEROQMRSQCcc18Ap/UG5Zz7rNTL1UCSmdV3zh2vwfCqYBdwegcUIiIi0S6a2wD8B/BZBQf/SWa20sxeNLNmNRmYpxCvwkKXAEQkXNoAFsKpTaV7jIuLIyMjg549ezJu3DiOHDlyxtGXHuFuwoQJFQ7PG4rhgGfNmsWjjz4KnD4ccOnRCMNh1KhRJQP31BZhSwDMbK6Z5QSYxgSxbg/gf4A7y1nkGeBcIAPYAfxvBduaaGbLzGxZXl5e1T9IufbgDQakBEBEwmVnjW+veCyAnJwcEhMTT+txrrCw8Iz2/MILL5CWllbu/FAkAFdddRX3338/cHoCEG6zZ88mOTm5xvYXCmFLAJxzlzrnegaY3qpoPTPrAMwEbnLOBey/0Tm30zlX6JwrAp4HBlQQx1TnXKZzLjMlJZT366sTIBGp24YMGcL69euZN28eF198Md///vdJT0+nsLCQyZMn079/f3r16sVzz3nNuZxzTJo0ibS0NK644gp27dpVsq3SZ+Dvvfceffv2pXfv3lxyySVs3ryZZ599lscff7zcroALCwvp3Lkzzjn27dtHvXr1yMrKOiXOadOmMWnSJBYuXMisWbOYPHkyGRkZJV0Bv/baawwYMICuXbuW290wwLRp0xg7diyXX345Xbp04Wc/+1nJvOnTp5Oenk7Pnj257777St7v1KkTu3fv5vDhw1xxxRX07t2bnj178uqrrwKwfPlyhg0bRr9+/Rg5ciQ7duw40z9LyERVPwBmlgy8AzzgnPu0guXaOueKS+8avEaFNaz4i62RAEWk7ikoKODdd9/l8ssvB2DJkiXk5OSQmprK1KlTadq0KUuXLuX48eMMHjyYESNG8Nlnn7F27VpWrVrFzp07SUtL47bbbjtlu3l5edxxxx1kZWWRmprK3r17ad68OXfddReNGjXi3nvvDRhPXFwcXbt2Zc2aNWzatIl+/frxySefMHDgQHJzcznvvPNYsGABABdeeCFXXXUVo0eP5tprrz3lMy1ZsoTZs2fz61//mrlz55b7+bOzs/nss8+oX78+3bp145577iEuLo777ruP5cuX06xZM0aMGMGbb77J1VdfXbLee++9R7t27XjnnXcAb6jk/Px87rnnHt566y1SUlJ49dVXefDBB3nxxRfP6G8TKpG6DfAaM8sFLgDeMbP3/VmTgPOA/y51i18rf50XSt0y+Jh/q+BK4GLgpzX9GU5WpakGQETqjqNHj5KRkUFmZiZnn302t99+OwADBgwgNTUVgDlz5vDyyy+TkZHBwIED2bNnD+vWrSMrK4vrr7+euLg42rVrx3e+c/qd3IsXL2bo0KEl22revHnQsQ0ZMoSsrCyysrJ44IEHWLBgAUuXLqV///5BrT927FgA+vXrx+bNmytc9pJLLqFp06YkJSWRlpbGli1bWLp0KcOHDyclJYX4+HhuuOGGklqIYunp6cydO5f77ruPTz75hKZNm7J27VpycnK47LLLyMjI4OGHHyY3Nzfozx0ukboLYCZeNX/Z9x8GHi5nnQmlnv8gfNEFS5cARKTuKW4DUFbDhg1LnjvnePLJJxk5cuQpy8yePfu0u7vKcs5Vukx5hgwZwrPPPsv27dt56KGHmDJlCvPmzWPo0KFBrV+/fn3Aq00oKCgIatnSywczeF7Xrl1Zvnw5s2fP5oEHHmDEiBFcc8019OjRg0WLFgUVZ02J5rsAotwuvPwpOcJxiIjUrJEjR/LMM8+Qn58PwFdffcXhw4cZOnQoM2bMoLCwkB07dvDxxx+ftu4FF1zA/Pnz2bRpEwB79+4FoHHjxhw8eLDC/Q4cOJCFCxdSr149kpKSyMjI4LnnnmPIkNNvxw5me1U1cOBA5s+fz+7duyksLGT69OmnDXm8fft2zjrrLG688UbuvfdeVqxYQbdu3cjLyytJAPLz81m9enVIYzsTSgDO2I+BRagIRSR8Qt3GKDTbmzBhAmlpafTt25eePXty5513UlBQwDXXXEOXLl1IT0/n7rvvPu3gCJCSksLUqVMZO3YsvXv35rrrrgPgyiuvZObMmeU2AgTvrLxjx44MGjQI8GoEDh48SHp6+mnLjh8/nilTptCnT5+SRoDV1bZtW373u99x8cUX07t3b/r27cuYMafe2LZq1SoGDBhARkYGjzzyCL/4xS9ITEzk9ddf57777qN3795kZGRU+46HULBgqjTqiszMTBfO+0BFRKoj0HjuIsEK9P0xs+XOucxAy+v0VUREJAZF1W2AIiIijzzyCK+99top740bNy7kQwS///77p9zLD5CamsrMmae1Ua+TdAlARCRK6BKAVIcuAYiI1GKxdFImoXMm3xslACIiUSIpKYk9e/YoCZAqcc6xZ88ekpKSqrSe2gCIiESJDh06kJubS2gHLpNYkJSURIcOHaq0jhIAEZEokZCQUNJFrki46RKAiIhIDFICICIiEoOUAIiIiMSgmOoHwMzygC0h3GRLYHcItxerVI6hoXIMDZVjaKgcQ6O65XiOcy4l0IyYSgBCzcyWldfBggRP5RgaKsfQUDmGhsoxNMJZjroEICIiEoOUAIiIiMQgJQDVMzXSAdQRKsfQUDmGhsoxNFSOoRG2clQbABERkRikGgAREZEYpAQgCGZ2uZmtNbP1ZnZ/gPlmZk/481eaWd9IxBntgijHG/zyW2lmC82sdyTijHaVlWOp5fqbWaGZXVuT8dUWwZSjmQ03s2wzW21m82s6xtogiP/rpmb2LzP73C/HWyMRZzQzsxfNbJeZ5ZQzPzzHGOecpgomIA7YAHQGEoHPgbQyy4wC3gUMGAT8O9JxR9sUZDleCDTzn39X5Xhm5VhquY+A2cC1kY472qYgv4/JwBrgbP91q0jHHW1TkOX4c+B//OcpwF4gMdKxR9MEDAX6AjnlzA/LMUY1AJUbAKx3zm10zp0AZgBjyiwzBnjZeRYDyWbWtqYDjXKVlqNzbqFz7lv/5WKgakNbxYZgvo8A9wBvALtqMrhaJJhy/D7wT+fcVgDnnMrydMGUowMam5kBjfASgIKaDTO6Oeey8MqlPGE5xigBqFx7YFup17n+e1VdJtZVtYxux8t45VSVlqOZtQeuAZ6twbhqm2C+j12BZmY2z8yWm9lNNRZd7RFMOT4FdAe2A6uAnzjnimomvDojLMcYDQdcOQvwXtlbJ4JZJtYFXUZmdjFeAnBRWCOqnYIpxz8C9znnCr2TLgkgmHKMB/oBlwANgEVmttg591W4g6tFginHkUA28B3gXOADM/vEOXcgzLHVJWE5xigBqFwu0LHU6w54mWxVl4l1QZWRmfUCXgC+65zbU0Ox1SbBlGMmMMM/+LcERplZgXPuzRqJsHYI9v96t3PuMHDYzLKA3oASgJOCKcdbgUeddzF7vZltAs4HltRMiHVCWI4xugRQuaVAFzNLNbNEYDwwq8wys4Cb/Jaag4D9zrkdNR1olKu0HM3sbOCfwA90llWuSsvROZfqnOvknOsEvA78UAf/0wTzf/0WMMTM4s3sLGAg8EUNxxntginHrXi1KJhZa6AbsLFGo6z9wnKMUQ1AJZxzBWY2CXgfr8Xri8651WZ2lz//WbyW1qOA9cARvIxXSgmyHH8JtACe9s9eC5wGEzlFkOUolQimHJ1zX5jZe8BKoAh4wTkX8DatWBXk9/E3wDQzW4VXlX2fc06jBJZiZtOB4UBLM8sF/h+QAOE9xqgnQBERkRikSwAiIiIxSAmAiIhIDFICICIiEoOUAIiIiMQgJQAiIiIxSAmAiIhIDFICICIiEoOUAIhI2JhZf3/88iQza+iPB98z0nGJiDoCEpEwM7OHgSS8AXVynXO/i3BIIoISABEJM7+P+KXAMeBC51xhhEMSEXQJQETCrznQCGiMVxMgIlFANQAiElZmNguYAaQCbZ1zkyIckoig0QBFJIzM7Ca8UR3/bmZxwEIz+45z7qNIxyYS61QDICIiEoPUBkBERCQGKQEQERGJQUoAREREYpASABERkRikBEBERCQGKQEQERGJQUoAREREYpASABERkRj0/wGXgqACHoBSNgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(8,5))\n",
    "S=plt.scatter(xs,ys_noise,s=5, marker='*',c='black')\n",
    "Predict_without_noise=plt.plot(x,y,c='red')\n",
    "Predict_with_noise=plt.plot(x,y_n,c='yellow')\n",
    "# Raw=plt.plot(x,F(x), c='g')\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "Samples=mpatches.Patch(color='black', label='Samples Points')\n",
    "Predict_without_noise=mpatches.Patch(color='red', label='Predict_without_noise')\n",
    "Predict_with_noise=mpatches.Patch(color='yellow', label='Predict_with_noise')\n",
    "plt.title('Successive_Linear_Interpolation with/withour noise')\n",
    "plt.legend(handles=[Samples,Predict_without_noise,Predict_with_noise],loc='best')\n",
    "# plt.savefig('T_R_100.jpg')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 946,
   "id": "451c9a1f",
   "metadata": {},
   "outputs": [],
   "source": [
    "LSI=Lagrange_Polynomial_Interpolation(x,y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 947,
   "id": "e316aeec",
   "metadata": {},
   "outputs": [],
   "source": [
    "for k in range(xs.shape[0]):\n",
    "    ys[k]=LSI.predict(xs[k])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 966,
   "id": "0e9e85eb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1c21a63d1c0>]"
      ]
     },
     "execution_count": 966,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAYKUlEQVR4nO3deXBddd3H8fe3C6ULpS1N00iRpBDowiI1Im0VkVJEYCgqbiO041b3bRwecHBUdJxBHsfBKo5WYZ44IogLUrEqNYgKIg8p3aWY2hYeaGnSolIs3b/PH78TkqZZbpJ77++cez6vmTPn3CW5Hw7pJye/+7vnmLsjIiLZMyR2ABERGRgVuIhIRqnARUQySgUuIpJRKnARkYwaVs4XmzhxotfW1pbzJUVEMm/lypU73b2q6/1lLfDa2lqam5vL+ZIiIplnZk91d7+GUEREMkoFLiKSUSpwEZGMUoGLiGSUClxEJKNU4CIiGaUCFxHJqGwU+K9/DTfdFDuFiEiqZKPAm5rgK18BnbtcRORl2Sjw2lrYswd27oydREQkNbJT4ABbt8ZMISKSKtkq8C1bosYQEUmTbBW4jsBFRF6WjQIfOxYmTFCBi4h0ko0Ch3AUriEUEZGXZafA6+p0BC4i0kl2Cry2NhS45oKLiABZK/C9e6G1NXYSEZFUyE6B19WFtcbBRUSAAgvczLaa2TozW21mzcl9E8xshZm1JOvxJU2qqYQiIkfozxH4G939Ve7ekNy+Hmhy93qgKbldOiefHNYqcBERYHBDKAuAxmS7Ebhy0Gl6M2YMTJyoAhcRSRRa4A7cb2YrzWxxcl+1u28HSNaTuvtCM1tsZs1m1tzW1ja4tHV1GgMXEUkMK/B5c919m5lNAlaY2cZCX8DdlwJLARoaGgY3B7C2FtasGdS3EBGpFAUdgbv7tmTdCtwDnAvsMLMagGRd+vl9tbXw1FNw+HDJX0pEJO36LHAzG21mx7VvAxcD64FlwKLkaYuAe0sV8mV1dbBvH+zYUfKXEhFJu0KGUKqBe8ys/fk/dvffmtljwN1m9n7gaeDtpYuZ6Hxa2Zqakr+ciEia9Vng7r4ZOLub+3cB80oRqked54LPmVPWlxYRSZvsfBITNBdcRKSTbBX4qFFQXa2phCIiZK3AQecFFxFJZK/Ap06FzZtjpxARiS57BX7qqWEu+P79sZOIiESVvQKvrw8f5NEbmSKSc9kr8FNPDeuWlrg5REQiy26Bb9oUN4eISGTZK/CJE+H441XgIpJ72Stws3AUrgIXkZzLXoFDKHCNgYtIzmW3wLduhQMHYicREYkmmwVeXw+HDoX54CIiOZXNAtdMFBGRjBe4xsFFJMeyWeCTJsFxx+kIXERyLZsFrqmEIiIZLXBQgYtI7mW7wDdvhoMHYycREYkiuwVeXx/K++mnYycREYkiuwWuqYQiknPZL3BNJRSRnMpugU+eDKNH6whcRHIruwVuFsbBn3wydhIRkSiyW+AA06bBxo2xU4iIRJHtAp8+PZyV8KWXYicRESm77Be4u4ZRRCSXsl/gAE88ETeHiEgEBRe4mQ01s1Vmdl9ye4KZrTCzlmQ9vnQxe1BfD0OGqMBFJJf6cwT+KaBzU14PNLl7PdCU3C6vESNg6lQVuIjkUkEFbmZTgMuAH3S6ewHQmGw3AlcWNVmhpk/XTBQRyaVCj8BvAf4LONzpvmp33w6QrCcVN1qBpk+Hv/9dJ7USkdzps8DN7HKg1d1XDuQFzGyxmTWbWXNbW9tAvkXvpk+H/fthy5bif28RkRQr5Ah8LnCFmW0F7gIuNLMfATvMrAYgWbd298XuvtTdG9y9oaqqqkixO5k2Law1Di4iOdNngbv759x9irvXAu8CHnD3q4FlwKLkaYuAe0uWsjeaSigiOTWYeeA3AfPNrAWYn9wuv+OPh5oaFbiI5M6w/jzZ3R8EHky2dwHzih9pAKZPV4GLSO5k+5OY7dqnErrHTiIiUjaVUeDTpsELL8D27bGTiIiUTWUUuN7IFJEcUoGLiGRUZRR4TQ2MGwcbNsROIiJSNpVR4GZw5pmwdm3sJCIiZVMZBQ5w1lmwbh0cPtz3c0VEKkBlFfju3fDUU7GTiIiURWUVOGgYRURyo3IK/IwzwloFLiI5UTkFPmYMnHKKClxEcqNyChzCMIoKXERyovIKvKUF9uyJnUREpOQqr8Dd4W9/i51ERKTkKq/AQcMoIpILlVXgU6fCqFEqcBHJhcoq8CFDwnRCFbiI5EBlFTh0zETRxR1EpMJVZoHv2qWLO4hIxavMAgdYsyZuDhGREqu8Aj/77LBetSpuDhGREqu8Ah83Dk49FVaujJ1ERKSkKq/AAV79amhujp1CRKSkKrPAGxrg6aehtTV2EhGRkqncAgcNo4hIRavMAp81K6xV4CJSwSqzwMeOhdNO0zi4iFS0yixwCMMoKnARqWCVXeDPPgvPPRc7iYhISfRZ4GZ2rJn9r5mtMbMNZnZjcv8EM1thZi3Jenzp4/aD3sgUkQpXyBH4PuBCdz8beBVwiZmdB1wPNLl7PdCU3E6Pc84BMw2jiEjF6rPAPXgxuTk8WRxYADQm9zcCV5Yi4ICNGQPTpqnARaRiFTQGbmZDzWw10AqscPdHgWp33w6QrCf18LWLzazZzJrb2tqKFLtADQ0aQhGRilVQgbv7IXd/FTAFONfMzij0Bdx9qbs3uHtDVVXVAGMOUENDOK3stm3lfV0RkTLo1ywUd/8X8CBwCbDDzGoAknX6Prf+mteE9V//GjeHiEgJFDILpcrMxiXbI4GLgI3AMmBR8rRFwL0lyjhws2bBMcfAI4/ETiIiUnTDCnhODdBoZkMJhX+3u99nZo8Ad5vZ+4GngbeXMOfAjBgRhlH+8pfYSUREiq7PAnf3tcA53dy/C5hXilBFNWcOLFkCe/fCscfGTiMiUjSV+0nMdnPnwv798PjjsZOIiBRV5Rf47NlhrWEUEakwlV/g1dVwyikqcBGpOJVf4BCGUR5+GNxjJxERKZp8FPicOeHyaps3x04iIlI0+Slw0DCKiFSUfBT4zJnhKj0qcBGpIPko8CFDwmyUhx+OnUREpGjyUeAQhlHWr4d//St2EhGRoshPgb/hDWEWyp//HDuJiEhR5KfAzzsvfJT+gQdiJxERKYr8FPiIEWE+uApcRCpEfgoc4MILYe1a2LkzdhIRkUHLX4EDPPhg1BgiIsWQrwJ/9avDxY41jCIiFSBfBT58OJx/vgpcRCpCvgocwjDKk0/Cs8/GTiIiMij5LHCAP/whbg4RkUHKX4GffTaMH68CF5HMy1+BDxkCF1wATU06P7iIZFr+Chxg/nx46qkwFi4iklH5LPBLLw3r5cvj5hARGYR8FvjJJ4dzhKvARSTD8lngEI7C//Qn2L07dhIRkQHJd4EfOAC//33sJCIiA5LfAp87N1xmTcMoIpJR+S3w4cPh4otDgWs6oYhkUH4LHMIwyrZt4RSzIiIZk+8Cf/Obw1rDKCKSQX0WuJmdZGZ/MLMnzGyDmX0quX+Cma0ws5ZkPb70cYts8uRwitn77oudRESk3wo5Aj8IfNbdpwPnAR8zsxnA9UCTu9cDTcnt7FmwAB55JAyliIhkSJ8F7u7b3f3xZHs38ARwIrAAaEye1ghcWaKMpfW2t4U3Me+5J3YSEZF+6dcYuJnVAucAjwLV7r4dQskDk3r4msVm1mxmzW1tbYOMWwIzZsD06fDzn8dOIiLSLwUXuJmNAX4OfNrdXyj069x9qbs3uHtDVVXVQDKW3lVXwR//CK2tsZOIiBSsoAI3s+GE8r7D3X+R3L3DzGqSx2uA7LbfVVfB4cPwy1/GTiIiUrBCZqEYcBvwhLt/o9NDy4BFyfYi4N7ixyuTM8+E+noNo4hIphRyBD4XuAa40MxWJ8ulwE3AfDNrAeYnt7PJLByFNzXBrl2x04iIFKSQWSgPubu5+1nu/qpkWe7uu9x9nrvXJ+vnyxG4ZN72Njh0CJYti51ERKQg+f4kZmezZkFdHfzkJ7GTiIgURAXezgze8x5YsUIf6hGRTFCBd3bNNWE2yh13xE4iItInFXhnp50Gs2dDY6NOMSsiqacC72rhQtiwAVatip1ERKRXKvCu3vlOGDEiHIWLiKSYCryr8ePhiivgxz8O18wUEUkpFXh3Fi6EnTvhN7+JnUREpEcq8O686U1QXQ3f/37sJCIiPVKBd2f4cPjgB+HXv4YtW2KnERHplgq8Jx/6EAwZAt/9buwkIiLdUoH3ZMqUcLm1226DvXtjpxEROYoKvDcf+1g4O6HOjyIiKaQC780b3xgut3brrbGTiIgcRQXeGzP46EfhscfCIiKSIirwvixcCGPHwte/HjuJiMgRVOB9GTsWPvIR+NnPoKUldhoRkZepwAvx6U+HueE33xw7iYjIy1TghZg8Gd73vnCCq2efjZ1GRARQgRfu2mvDxR6+8Y3YSUREABV44erq4F3vgu99T1euF5FUUIH3x+c+B3v2aCxcRFJBBd4fM2fC1VfDkiXwzDOx04hIzqnA++vGG+HQIfjyl2MnEZGcU4H3V11dmBd+++3w5JOx04hIjqnAB+KGG2DkSPj852MnEZEcU4EPxKRJ8NnPhk9nPvRQ7DQiklMq8IG69lo46aRwytmDB2OnEZEc6rPAzex2M2s1s/Wd7ptgZivMrCVZjy9tzBQaPRpuuQXWrtXpZkUkikKOwP8HuKTLfdcDTe5eDzQlt/PnLW8JF0D+whdg+/bYaUQkZ/oscHf/E/B8l7sXAI3JdiNwZXFjZYQZfOtb4ZJr114bO42I5MxAx8Cr3X07QLKeVLxIGVNfD9ddB3fcEa5iLyJSJiV/E9PMFptZs5k1t7W1lfrl4rjhBjjrLPjAB3SeFBEpm4EW+A4zqwFI1q09PdHdl7p7g7s3VFVVDfDlUm7ECPjhD0N5f/SjsdOISE4MtMCXAYuS7UXAvcWJk2Fnnw1f/CLcfTfcdVfsNCKSA4VMI7wTeAQ43cyeMbP3AzcB882sBZif3JbrroPXvhY+/GHYvDl2GhGpcMP6eoK7v7uHh+YVOUv2DRsGd94Js2bBVVfBww+Hj9yLiJSAPolZbHV18KMfwapV8IlPxE4jIhVMBV4Kl10WZqbcdhv84Aex04hIhVKBl8qNN8L8+eHUsw88EDuNiFQgFXipDB0aZqScfjq89a2wfn3fXyMi0g8q8FIaNw6WL4dRo+DSS2HbttiJRKSCqMBL7ZWvDB+x/+c/w5BKa4+feRIR6RcVeDmccw786lewZQvMmweVekoBESkrFXi5XHBBKPFNm+Cii2DnztiJRCTjVODlNG8eLFsWLob8utfB1q2xE4lIhqnAy23+fFixAnbsgNmzYfXq2IlEJKNU4DG8/vXhYsjDhsH554eZKiIi/aQCj2XmTHjkETjlFLj8cvjKV+Dw4dipRCRDVOAxTZkSTnh19dXhuppXXqkLQohIwVTgsY0aBY2NsGQJ/Pa3cOaZcP/9sVOJSAaowNPALJy58NFHw6c33/SmcHv37tjJRCTFVOBpcs45sHIlfPKT8O1vw4wZcM894B47mYikkAo8bUaOhG9+E/7yF5gwIZwI6/LLYcOG2MlEpD/27YPHHoNbb4VFi6Clpegv0ecVeSSS2bPD0fiSJeHUtGedBe99L3zpS+HNTxFJjz17YN26cCGX9mXNGti/Pzw+aRJccw3U1xf1Zc3L+Od5Q0ODNzc3l+31KsauXfDVr4ZhFTNYuBCuvRZOOy12MpF8cYdnnoEnnoC1a0NRr14NGzd2TAMeNy4Mh77mNWE591w46aTwb3eAzGyluzccdb8KPEO2boWbb4bbbw+/2a+6Cj7zGTjvvEH9cIhIF//5Tzj53D/+Ecr6b38L640b4cUXO543ZUoo687LK19Z9H+PKvBKsmMH3HILfOc78MILcMYZsHhxmE8+fnzsdCLpdvhwOCPotm2wfTs8+2w4ONqyBTZvDuuup30+8USYPj0sM2aE9cyZMHFiWSKrwCvR7t1w112wdCk0N8Oxx4brcb7jHWE9enTshCKl5R7Gn//5T3j++TDc+PzzRy+traGst22D556DQ4eO/D5Dh8LJJ4eLkk+dGtbt29Omwdixcf77EirwSrdqVbiI8s9+Fo7QR44MVwG67LIwr/wVr4idUCrB4cNw4EBY9u8/cj2Y+7p77KWXwlDGiy8evW7f/s9/ep9me8wxcMIJUFUV/g3U1By57rw9LL1zOlTgeXHoEPz5z/DTn8IvfhGONiB8wvPii2HuXJgzB6qr4+aUvrmHqWgvvdT/Ze/ejjLsXJLd3S70sQMHjj5yLTazULrDh4e/KMeMCcvo0b2vJ0w4cjnhhLAeObIi3h9SgeeRe3in/He/C8tDD3VMazrllFDk554byv3MM8MPvAzO/v1haOuFF8K683Zv93V97MUXQwkP9N/nkCEwYkRHGR5zzNHbhT7W27rY9w0dWtz/HxVCBS7haO7xx8OHhB5+OKx37Oh4/MQTw3zz00/vGP+bOhVqa8M5WyrVvn2FF21f5btvX2GvOXJkGFc97riOdeftMWPCcwa6DB9e2n0mZaUCl6O5hzd11q0Ly9q1Yb1pUxhb7Gzy5DBWWF0dlsmTw7qqCo4/PhRP12XEiOJlPXgwlOPevR1L+1DB3r3hiHX37o51T9vd3df+V0lfRo06unS7K9++HhszJtXjrZI+PRW4foryzCwcdZ94IlxyScf97mGa1ebNHdOqtmwJ4+nPPQfr14cj9wMH+v7+Pf3p3H6EePhwz8vBgx0FffBg///7hg3rKND2o9rjjgu/fNq3x4zp+IXTV+nqz3tJGRW4HM0sfPR30qTwIaHuuIepWzt3hiGErsu//33kG2ldZxgcOBBeZ8iQnpehQ8NwwLHH9r6MGNFRyJ2Luph/AYik0KAK3MwuAb4JDAV+4O43FSWVpJ9Zxzv+IhLFgM9GaGZDgVuBNwMzgHeb2YxiBRMRkd4N5nSy5wKb3H2zu+8H7gIWFCeWiIj0ZTAFfiLwf51uP5PcdwQzW2xmzWbW3NbWNoiXExGRzgZT4N19vOmoOYnuvtTdG9y9oaqqahAvJyIinQ2mwJ8BTup0ewqwbXBxRESkUIMp8MeAejOrM7NjgHcBy4oTS0RE+jLgaYTuftDMPg78jjCN8HZ314UbRUTKZFDzwN19ObC8SFlERKQfynouFDNrA54a4JdPBHYWMU65KX88Wc4Oyh9TWrKf7O5HzQIpa4EPhpk1d3cyl6xQ/niynB2UP6a0Zx/Mm5giIhKRClxEJKOyVOBLYwcYJOWPJ8vZQfljSnX2zIyBi4jIkbJ0BC4iIp2owEVEMiq1BW5m/21mG81srZndY2bjenjeJWb2pJltMrPryxyzR2b2djPbYGaHzazHaUhmttXM1pnZajNLzQVD+5E/dfvfzCaY2Qoza0nW43t4Xqr2fV/70oIlyeNrzWxWjJzdKSD7BWb272RfrzazL8TI2RMzu93MWs1sfQ+Pp3Pfu3sqF+BiYFiy/TXga908ZyjwD2AqcAywBpgRO3uSbTpwOvAg0NDL87YCE2PnHUj+tO5/4Gbg+mT7+u5+dtK27wvZl8ClwG8IZwI9D3g0du5+ZL8AuC921l7+G84HZgHre3g8lfs+tUfg7n6/u7dfyfavhLMddpXai0q4+xPu/mTsHANVYP607v8FQGOy3QhcGS9KwQrZlwuAH3rwV2CcmdWUO2g30vpzUDB3/xPwfC9PSeW+T22Bd/E+wm+/rgq6qETKOXC/ma00s8Wxw/RTWvd/tbtvB0jWk3p4Xpr2fSH7Mq37u9Bcs81sjZn9xsxmlida0aRy30e9Kr2Z/R6Y3M1DN7j7vclzbgAOAnd09y26ua9s8yILyV+Aue6+zcwmASvMbGNyNFByRcgfbf/3lr0f3ybavu9GIfsy6s97LwrJ9TjhfB4vmtmlwC+B+lIHK6JU7vuoBe7uF/X2uJktAi4H5nkyENVF1ItK9JW/wO+xLVm3mtk9hD9Hy1IiRcgfbf/3lt3MdphZjbtvT/7Mbe3he0Tb990oZF+m9SIqfeZy9xc6bS83s++Y2UR3T8OJogqRyn2f2iEUM7sEuA64wt339PC0TF9UwsxGm9lx7duEN267fRc8pdK6/5cBi5LtRcBRf02kcN8Xsi+XAQuTGRHnAf9uHyqKrM/sZjbZzCzZPpfQPbvKnnTg0rnvY7+L2tMCbCKMOa1Olu8m978CWN7l3eG/E94FvyF27k653kL4rb0P2AH8rmt+wrv2a5JlQ9byp3X/AycATUBLsp6QhX3f3b4EPgx8ONk24Nbk8XX0Mrsphdk/nuznNYRJCXNiZ+6S/05gO3Ag+bl/fxb2vT5KLyKSUakdQhERkd6pwEVEMkoFLiKSUSpwEZGMUoGLiGSUClxEJKNU4CIiGfX/l4tkp4pVzSAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(xs,ys,c='r')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4d375ad4",
   "metadata": {},
   "source": [
    "## Newton插值"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "058ad59b",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.set_printoptions(precision=5, suppress=True, linewidth=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "87a141e1",
   "metadata": {},
   "outputs": [],
   "source": [
    "def cal_dc(x,y, eps=1e-3):\n",
    "    Difference_coefficient=lambda a,b,x1,x2:(a-b)/(x[x1]-x[x2])\n",
    "    N=x.shape[0]\n",
    "    dc_matrix=np.zeros((N,N))\n",
    "    base_order=0\n",
    "    for i in range(N):\n",
    "        for j in range(i+1):\n",
    "            if j==0:\n",
    "                dc_matrix[i][j]=y[i]\n",
    "            else:\n",
    "                dc_matrix[i][j]=Difference_coefficient(dc_matrix[i][j-1],dc_matrix[i-1][j-1],i,i-j)\n",
    "                if abs(dc_matrix[i][j]-dc_matrix[i-1][j])<=eps:\n",
    "                    base_order=i-1\n",
    "                \n",
    "    for i in range(dc_matrix.shape[0]):\n",
    "            print('{}:'.format(y[i]), end='\\t')\n",
    "            for j in range(dc_matrix.shape[1]):\n",
    "                print(dc_matrix[i][j], end='\\t')\n",
    "            print('\\n')\n",
    "    print('base_order={}'.format(base_order))\n",
    "    return dc_matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "id": "d0d1489b",
   "metadata": {},
   "outputs": [],
   "source": [
    "x=np.array([0.4,0.55,0.65,0.80,0.90,1.05]).astype(np.float64)\n",
    "y=np.array([0.41075,0.57815,0.69675,0.88811,1.02652,1.25382]).astype(np.float64)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "3817a918",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.41075:\t0.41075\t0.0\t0.0\t0.0\t0.0\t0.0\t\n",
      "\n",
      "0.57815:\t0.57815\t1.116\t0.0\t0.0\t0.0\t0.0\t\n",
      "\n",
      "0.69675:\t0.69675\t1.1859999999999995\t0.2799999999999976\t0.0\t0.0\t0.0\t\n",
      "\n",
      "0.88811:\t0.88811\t1.275733333333333\t0.35893333333333377\t0.19733333333334047\t0.0\t0.0\t\n",
      "\n",
      "1.02652:\t1.02652\t1.3841000000000017\t0.4334666666666749\t0.21295238095240318\t0.03123809523812543\t0.0\t\n",
      "\n",
      "1.25382:\t1.25382\t1.515333333333332\t0.5249333333333217\t0.22866666666661706\t0.031428571428427754\t0.0002930402927728068\t\n",
      "\n",
      "base_order=4\n"
     ]
    }
   ],
   "source": [
    "dc=cal_dc(x,y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "590a45cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "def cal_y(dc,x,y,target,order):\n",
    "    error=abs(dc[dc.shape[0]-1][dc.shape[0]-1]*reduce(lambda a,b: a*b, [target-x[k] for k in range(dc.shape[0])]))\n",
    "    print('error:{}'.format(error))\n",
    "    return sum([dc[i][i]*(1 if i==0 else reduce(lambda a,b: a*b, [target-x[k] for k in range(i)]) ) for i in range(order)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "4766f406",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "error:4.01693316910874e-09\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.6319144055039999"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cal_y(dc,x,y,0.596,4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "id": "bcdc0f78",
   "metadata": {},
   "outputs": [],
   "source": [
    "class Newton_Linear_Interpolation:\n",
    "    \"\"\"Newton_Linear_Interpolation 逐次线性插值\n",
    "\n",
    "    Author: Junno\n",
    "    Date: 2022-03-03\n",
    "    \"\"\"\n",
    "    def __init__(self,x,y,eps=1e-3,check=False):\n",
    "        \"\"\"__init__ \n",
    "\n",
    "        Args:\n",
    "            x ([np.array]): [x of points]\n",
    "            y ([np.array]): [y of points]\n",
    "            eps ([float], optional): [threshold to stop process]. Defaults to 1e-3.\n",
    "            check ([bool], optional): [whether to show dc_matrix]. Defaults to False.\n",
    "        \"\"\"\n",
    "        self.N=x.shape[0]\n",
    "        self.dc_matrix=np.zeros((self.N,self.N))\n",
    "        self.base_order=0\n",
    "        self.eps=eps\n",
    "        self.x=x\n",
    "        self.y=y\n",
    "        \n",
    "        self.Difference_coefficient=lambda a,b,x1,x2:(a-b)/(self.x[x1]-self.x[x2])\n",
    "        \n",
    "        for i in range(self.N):\n",
    "            for j in range(i+1):\n",
    "                if j==0:\n",
    "                    self.dc_matrix[i][j]=y[i]\n",
    "                else:\n",
    "                    self.dc_matrix[i][j]=self.Difference_coefficient(self.dc_matrix[i][j-1],self.dc_matrix[i-1][j-1],i,i-j)\n",
    "                    if abs(self.dc_matrix[i][j]-self.dc_matrix[i-1][j])<=self.eps:\n",
    "                        self.base_order=i-1\n",
    "\n",
    "        if self.base_order == 0:\n",
    "            self.base_order=self.N-1\n",
    "\n",
    "        if check:    \n",
    "            self.show_dc()\n",
    "    \n",
    "    def calulate(self,x,order,show_error=True):\n",
    "        \"\"\"__init__ \n",
    "\n",
    "        Args:\n",
    "            x ([float]): [x to calculate y]\n",
    "            order ([itn]): [demanding order]\n",
    "            show_error ([bool], optional): [whether to error]. Defaults to True.\n",
    "            \n",
    "        Return:\n",
    "            y ([float]): [predict y value]\n",
    "        \"\"\"\n",
    "        assert order<self.N,'order should be less than the number of points'\n",
    "        t_order = min(order,self.base_order+1)\n",
    "        error=abs(self.dc_matrix[t_order][t_order]*reduce(lambda a,b: a*b, [x-self.x[k] for k in range(self.dc_matrix.shape[0])]))\n",
    "        if show_error:\n",
    "            print('error:{}'.format(error))\n",
    "        return sum([self.dc_matrix[i][i]*(1 if i==0 else reduce(lambda a,b: a*b, [x-self.x[k] for k in range(i)]) ) for i in range(t_order)])\n",
    "    \n",
    "    def add_point(self, x, y):\n",
    "        \"\"\"add_point \n",
    "\n",
    "        Args:\n",
    "            x ([np.array]): [x of points]\n",
    "            y ([np.array]): [y of points]\n",
    "        \"\"\"\n",
    "        M = x.shape[0]\n",
    "        self.x = np.concatenate((self.x, x), axis=0)\n",
    "        self.y = np.concatenate((self.y, y), axis=0)\n",
    "        self.N = self.N+M\n",
    "        self.dc_matrix=np.pad(self.dc_matrix,((0,M),(0,M)))\n",
    "        for i in range(self.base_order+1,self.N):\n",
    "            for j in range(i+1):\n",
    "                if j==0:\n",
    "                    self.dc_matrix[i][j]=self.y[i]\n",
    "                else:\n",
    "                    self.dc_matrix[i][j]=self.Difference_coefficient(self.dc_matrix[i][j-1],self.dc_matrix[i-1][j-1],i,i-j)\n",
    "                    if abs(self.dc_matrix[i][j]-self.dc_matrix[i-1][j])<=self.eps:\n",
    "                        self.base_order=i-1\n",
    "                        \n",
    "        self.show_dc()\n",
    "        \n",
    "    def show_dc(self):\n",
    "        \"\"\"show difference coefficient matrix\n",
    "        \"\"\"\n",
    "        for i in range(self.dc_matrix.shape[0]):\n",
    "            print('{}:'.format(self.y[i]), end='\\t')\n",
    "            for j in range(self.dc_matrix.shape[1]):\n",
    "                print(self.dc_matrix[i][j], end='\\t')\n",
    "            print('\\n')\n",
    "        print('base_order={}'.format(self.base_order))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "id": "71ec9a50",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.41075:\t0.41075\t0.0\t0.0\t0.0\t0.0\t0.0\t\n",
      "\n",
      "0.57815:\t0.57815\t1.116\t0.0\t0.0\t0.0\t0.0\t\n",
      "\n",
      "0.69675:\t0.69675\t1.1859999999999995\t0.2799999999999976\t0.0\t0.0\t0.0\t\n",
      "\n",
      "0.88811:\t0.88811\t1.275733333333333\t0.35893333333333377\t0.19733333333334047\t0.0\t0.0\t\n",
      "\n",
      "1.02652:\t1.02652\t1.3841000000000017\t0.4334666666666749\t0.21295238095240318\t0.03123809523812543\t0.0\t\n",
      "\n",
      "1.25382:\t1.25382\t1.515333333333332\t0.5249333333333217\t0.22866666666661706\t0.031428571428427754\t0.0002930402927728068\t\n",
      "\n",
      "base_order=4\n"
     ]
    }
   ],
   "source": [
    "NSL=Newton_Linear_Interpolation(x,y,check=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "id": "04088d73",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "error:4.01693316910874e-09\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.6319175080796159"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "NSL.calulate(0.596,5,show_error=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "75b02ba5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.41075:\t0.41075\t0.0\t0.0\t0.0\t0.0\t0.0\t\n",
      "\n",
      "0.57815:\t0.57815\t1.116\t0.0\t0.0\t0.0\t0.0\t\n",
      "\n",
      "0.69675:\t0.69675\t1.1859999999999995\t0.2799999999999976\t0.0\t0.0\t0.0\t\n",
      "\n",
      "0.88811:\t0.88811\t1.275733333333333\t0.35893333333333377\t0.19733333333334047\t0.0\t0.0\t\n",
      "\n",
      "1.02652:\t1.02652\t1.3841000000000017\t0.4334666666666749\t0.21295238095240318\t0.03123809523812543\t0.0\t\n",
      "\n",
      "1.25382:\t1.25382\t1.515333333333332\t0.5249333333333217\t0.22866666666661706\t0.031428571428427754\t0.0002930402927728068\t\n",
      "\n",
      "base_order=4\n"
     ]
    }
   ],
   "source": [
    "NSL.show_dc()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "f52452af",
   "metadata": {},
   "outputs": [],
   "source": [
    "x=np.array([0.4,0.55,0.65,0.80]).astype(np.float64)\n",
    "y=np.array([0.41075,0.57815,0.69675,0.88811]).astype(np.float64)\n",
    "\n",
    "ax=np.array([0.90,1.05]).astype(np.float64)\n",
    "ay=np.array([1.02652,1.25382]).astype(np.float64)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "2ad66bcb",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.41075:\t0.41075\t0.0\t0.0\t0.0\t\n",
      "\n",
      "0.57815:\t0.57815\t1.116\t0.0\t0.0\t\n",
      "\n",
      "0.69675:\t0.69675\t1.1859999999999995\t0.2799999999999976\t0.0\t\n",
      "\n",
      "0.88811:\t0.88811\t1.275733333333333\t0.35893333333333377\t0.19733333333334047\t\n",
      "\n",
      "base_order=3\n"
     ]
    }
   ],
   "source": [
    "NSL=Newton_Linear_Interpolation(x,y,check=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "2be9321e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.41075:\t0.41075\t0.0\t0.0\t0.0\t0.0\t0.0\t\n",
      "\n",
      "0.57815:\t0.57815\t1.116\t0.0\t0.0\t0.0\t0.0\t\n",
      "\n",
      "0.69675:\t0.69675\t1.1859999999999995\t0.2799999999999976\t0.0\t0.0\t0.0\t\n",
      "\n",
      "0.88811:\t0.88811\t1.275733333333333\t0.35893333333333377\t0.19733333333334047\t0.0\t0.0\t\n",
      "\n",
      "1.02652:\t1.02652\t1.3841000000000017\t0.4334666666666749\t0.21295238095240318\t0.03123809523812543\t0.0\t\n",
      "\n",
      "1.25382:\t1.25382\t1.515333333333332\t0.5249333333333217\t0.22866666666661706\t0.031428571428427754\t0.0002930402927728068\t\n",
      "\n",
      "base_order=4\n"
     ]
    }
   ],
   "source": [
    "NSL.add_point(ax,ay)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "9ec90bac",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "error:3.838179646586846e-06\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.629486"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "NSL.calulate(0.596,2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "554af46a",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3c3546bf",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "ca97c54a",
   "metadata": {},
   "source": [
    "## 差分Newton插值多项式"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "ea1131a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 二项式系数\n",
    "Binomial_coefficient=lambda j,k:reduce(lambda x,y:x*y,[j-i for i in range(k)])/reduce(lambda x,y:x*y, [i for i in range(1,k+1)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "7c2c4239",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.125"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Binomial_coefficient(0.5,2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "id": "37d2d157",
   "metadata": {},
   "outputs": [],
   "source": [
    "x=np.array([0.4,0.55,0.65,0.80]).astype(np.float64)\n",
    "y=np.array([0.41075,0.57815,0.69675,0.88811]).astype(np.float64)\n",
    "temp=np.zeros((y.shape[0],y.shape[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16bfdf2d",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "id": "fd0a2ca5",
   "metadata": {},
   "outputs": [],
   "source": [
    "temp[:,0]=y\n",
    "for j in range(1,y.shape[0]):\n",
    "    for i in range(y.shape[0]-j):\n",
    "        temp[i][j]=temp[i+1][j-1]-temp[i][j-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "id": "e2e0069d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.41075,  0.1674 , -0.0488 ,  0.12156],\n",
       "       [ 0.57815,  0.1186 ,  0.07276,  0.     ],\n",
       "       [ 0.69675,  0.19136,  0.     ,  0.     ],\n",
       "       [ 0.88811,  0.     ,  0.     ,  0.     ]])"
      ]
     },
     "execution_count": 86,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "temp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "id": "0edf0131",
   "metadata": {},
   "outputs": [],
   "source": [
    "y=np.array([0.41075,0.57815,0.6532,0.69675,0.88811]).astype(np.float64)\n",
    "temp=np.zeros((y.shape[0],y.shape[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "id": "6d7ab71b",
   "metadata": {},
   "outputs": [],
   "source": [
    "temp[:,0]=y\n",
    "for j in range(1,y.shape[0]):\n",
    "    for i in range(y.shape[0]-j):\n",
    "        temp[i][j]=temp[i+1][j-1]-temp[i][j-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "id": "3b5bc799",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.41075,  0.1674 , -0.09235,  0.06085,  0.11846],\n",
       "       [ 0.57815,  0.07505, -0.0315 ,  0.17931,  0.     ],\n",
       "       [ 0.6532 ,  0.04355,  0.14781,  0.     ,  0.     ],\n",
       "       [ 0.69675,  0.19136,  0.     ,  0.     ,  0.     ],\n",
       "       [ 0.88811,  0.     ,  0.     ,  0.     ,  0.     ]])"
      ]
     },
     "execution_count": 97,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "temp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b31fccc6",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 391,
   "id": "100132cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "class DNewton_Linear_Interpolation:\n",
    "    \"\"\"DNewton_Linear_Interpolation 逐次线性插值\n",
    "\n",
    "    Author: Junno\n",
    "    Date: 2022-03-03\n",
    "    \"\"\"\n",
    "    def __init__(self,x,y,h,eps=1e-3,check=False):\n",
    "        \"\"\"__init__ \n",
    "\n",
    "        Args:\n",
    "            x ([np.array]): [x of points]\n",
    "            y ([np.array]): [y of points]\n",
    "            eps ([float], optional): [threshold to stop process]. Defaults to 1e-3.\n",
    "            check ([bool], optional): [whether to show dc_matrix]. Defaults to False.\n",
    "        \"\"\"\n",
    "        self.N=y.shape[0]\n",
    "        self.base_order=0\n",
    "        self.eps=eps\n",
    "        self.x=x\n",
    "        self.y=y\n",
    "        self.h=h\n",
    "        \n",
    "        self.binomial_coefficient=lambda j,k:reduce(lambda x,y:x*y,[j-i for i in range(k)])/reduce(lambda x,y:x*y, [i for i in range(1,k+1)]) if k>0 else 1\n",
    "        \n",
    "        self.forward_delta_matrix=self.cal_forward_(self.y,self.N)\n",
    "                    \n",
    "        self.backward_nabla_matrix=self.cal_backward_(self.y,self.N)\n",
    "\n",
    "        if self.base_order == 0:\n",
    "            self.base_order=self.N-1\n",
    "    \n",
    "    def cal_forward_(self,y,N):\n",
    "        temp=np.zeros((N,N))\n",
    "        temp[:,0]=y\n",
    "        for j in range(1,N):\n",
    "            for i in range(N-j):\n",
    "                temp[i][j]=temp[i+1][j-1]-temp[i][j-1]\n",
    "        return temp \n",
    "        \n",
    "    def cal_backward_(self,y,N):\n",
    "        temp=np.zeros((N,N))\n",
    "        temp[:,0]=y\n",
    "        for j in range(1,N):\n",
    "            for i in range(j,N):\n",
    "                temp[i][j]=temp[i-1][j-1]-temp[i][j-1]\n",
    "        return temp \n",
    "        \n",
    "    def calulate(self,x,eps=1e-3,direction=None,show_error=True):\n",
    "        \"\"\"__init__ \n",
    "\n",
    "        Args:\n",
    "            x ([float]): [x to calculate y]\n",
    "            order ([itn]): [demanding order]\n",
    "            show_error ([bool], optional): [whether to error]. Defaults to True.\n",
    "            \n",
    "        Return:\n",
    "            y ([float]): [predict y value]\n",
    "        \"\"\"\n",
    "        if direction==None:\n",
    "        # half to be forward and half to be backward\n",
    "            direction='forward' if x<self.x[0]+self.N/2*self.h else 'backward'\n",
    "        # find nearest points index\n",
    "        ind=np.argmin(abs(self.x-x))\n",
    "        \n",
    "        if direction == 'forward':\n",
    "            t=(x-self.x[ind])/self.h\n",
    "            \n",
    "            count=np.sum(np.abs(self.forward_delta_matrix[ind])>eps)\n",
    "            ans= sum([self.binomial_coefficient(t,k)*self.forward_delta_matrix[ind][k] for k in range(count)])\n",
    "\n",
    "            if show_error:\n",
    "                ny=np.sort(np.concatenate((self.y,[ans]),axis=0), axis=0)\n",
    "                nm=self.cal_forward_(ny,self.N+1)\n",
    "                error=np.abs(self.binomial_coefficient(t,count)*np.max(np.abs(nm[:,count])))\n",
    "                print(\"MAX INTERPOLATION ERROR:{}\".format(error))\n",
    "            return ans\n",
    "            \n",
    "        else:\n",
    "            t=-(x-self.x[ind])/self.h       \n",
    "            count=np.sum(np.abs(self.backward_nabla_matrix[ind])>eps)\n",
    "            ans=sum([(-1)**k*self.binomial_coefficient(-t,k)*self.backward_nabla_matrix[ind][k] for k in range(count)])\n",
    "            if show_error:\n",
    "                ny=np.sort(np.concatenate((self.y,[ans]),axis=0), axis=0)\n",
    "                nm=self.cal_backward_(ny,self.N+1)\n",
    "                error=np.abs(self.binomial_coefficient(-t,count)*np.max(np.abs(nm[:,count])))\n",
    "                print(\"MAX INTERPOLATION ERROR:{}\".format(error))\n",
    "                \n",
    "            return ans\n",
    "    \n",
    "    def add_point(self, x, y):\n",
    "        \"\"\"add_point \n",
    "\n",
    "        Args:\n",
    "            x ([np.array]): [x of points]\n",
    "            y ([np.array]): [y of points]\n",
    "        \"\"\"\n",
    "        M = x.shape[0]\n",
    "        self.x = np.concatenate((self.x, x), axis=0)\n",
    "        self.y = np.concatenate((self.y, y), axis=0)\n",
    "        self.N = self.N+M\n",
    "                "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 416,
   "id": "66edfc8a",
   "metadata": {},
   "outputs": [],
   "source": [
    "x=np.arange(0,2,0.2)\n",
    "F=lambda x:np.sin(x)+np.cos(x**2)\n",
    "y=F(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 417,
   "id": "b607161f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([0. , 0.2, 0.4, 0.6, 0.8, 1. , 1.2, 1.4, 1.6, 1.8]),\n",
       " array([ 1.     ,  1.19787,  1.37665,  1.50054,  1.51945,  1.38177,  1.06246,  0.606  ,  0.16398,\n",
       "        -0.02131]))"
      ]
     },
     "execution_count": 417,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x,y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 418,
   "id": "bad26427",
   "metadata": {},
   "outputs": [],
   "source": [
    "DN=DNewton_Linear_Interpolation(x,y,0.2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 303,
   "id": "b5b51bff",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.     ,  0.19867, -0.00792, -0.0076 ,  0.00062,  0.00028, -0.00004, -0.00001,  0.     ,\n",
       "         0.     ],\n",
       "       [ 0.19867,  0.19075, -0.01552, -0.00699,  0.0009 ,  0.00024, -0.00005, -0.00001,  0.     ,\n",
       "         0.     ],\n",
       "       [ 0.38942,  0.17522, -0.02251, -0.00609,  0.00114,  0.0002 , -0.00005, -0.00001,  0.     ,\n",
       "         0.     ],\n",
       "       [ 0.56464,  0.15271, -0.0286 , -0.00495,  0.00134,  0.00014, -0.00006,  0.     ,  0.     ,\n",
       "         0.     ],\n",
       "       [ 0.71736,  0.12411, -0.03355, -0.00361,  0.00148,  0.00008,  0.     ,  0.     ,  0.     ,\n",
       "         0.     ],\n",
       "       [ 0.84147,  0.09057, -0.03716, -0.00213,  0.00157,  0.     ,  0.     ,  0.     ,  0.     ,\n",
       "         0.     ],\n",
       "       [ 0.93204,  0.05341, -0.03929, -0.00056,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,\n",
       "         0.     ],\n",
       "       [ 0.98545,  0.01412, -0.03985,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,\n",
       "         0.     ],\n",
       "       [ 0.99957, -0.02573,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,\n",
       "         0.     ],\n",
       "       [ 0.97385,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,\n",
       "         0.     ]])"
      ]
     },
     "execution_count": 303,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "DN.forward_delta_matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 304,
   "id": "67aa7804",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,\n",
       "         0.     ],\n",
       "       [ 0.19867, -0.19867,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,\n",
       "         0.     ],\n",
       "       [ 0.38942, -0.19075, -0.00792,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,\n",
       "         0.     ],\n",
       "       [ 0.56464, -0.17522, -0.01552,  0.0076 ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,\n",
       "         0.     ],\n",
       "       [ 0.71736, -0.15271, -0.02251,  0.00699,  0.00062,  0.     ,  0.     ,  0.     ,  0.     ,\n",
       "         0.     ],\n",
       "       [ 0.84147, -0.12411, -0.0286 ,  0.00609,  0.0009 , -0.00028,  0.     ,  0.     ,  0.     ,\n",
       "         0.     ],\n",
       "       [ 0.93204, -0.09057, -0.03355,  0.00495,  0.00114, -0.00024, -0.00004,  0.     ,  0.     ,\n",
       "         0.     ],\n",
       "       [ 0.98545, -0.05341, -0.03716,  0.00361,  0.00134, -0.0002 , -0.00005,  0.00001,  0.     ,\n",
       "         0.     ],\n",
       "       [ 0.99957, -0.01412, -0.03929,  0.00213,  0.00148, -0.00014, -0.00005,  0.00001,  0.     ,\n",
       "         0.     ],\n",
       "       [ 0.97385,  0.02573, -0.03985,  0.00056,  0.00157, -0.00008, -0.00006,  0.00001,  0.     ,\n",
       "        -0.     ]])"
      ]
     },
     "execution_count": 304,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "DN.backward_nabla_matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 427,
   "id": "5cfa2efe",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.  0.2 0.4 0.6 0.8 1.  1.2 1.4 1.6 1.8]\n",
      "MAX INTERPOLATION ERROR:0.05321699743238864\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.18476646350550527"
      ]
     },
     "execution_count": 427,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d1='forward'\n",
    "d2='backward'\n",
    "v=1.59\n",
    "d=d1\n",
    "print(x)\n",
    "DN.calulate(v,show_error=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 428,
   "id": "a48e61f5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.0025930434075897846"
      ]
     },
     "execution_count": 428,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F(v)-DN.calulate(v,show_error=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 413,
   "id": "eb8f4ebc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.11436409003036918"
      ]
     },
     "execution_count": 413,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F(v)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 233,
   "id": "8b4e6a88",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2"
      ]
     },
     "execution_count": 233,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sum(x>0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45ce5ba7",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
